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Abstract 


The present research is concerned with the determination of the local Nusselt 
number for fluid flow and heat transfer from flat and ribbed surfaces. The deter- 
mination of local Nusselt number for a surface with known temperatures requires 
the temperature in the vicinity of the wail. This is equivalent to solving the 
energy equation with proper initial and boundary conditions. In turn, it calls 
for the accurate measurement of surface temperatures or surface heat fluxes. In 
practice, the surface temperature or surface heat flux measurements are generally 
prone to substantial amount of spatial and temporal errors. The present work 
shows that the local Nusselt number can be estimated with a good degree of ac- 
curacy, even from the noisy surface temperature data, if additional high quality 
data is available within a favorable region of the flow field. The improvement has 
been effected by an inverse forced convection formulation in the present work and 
solved by the conjugate-gradient approach. 

The present investigation is restricted to steady laminar flow of an incompress- 
ible Newtonian fluid with negligible viscous dissipation and mechanical work. The 
thermal field is allowed to evolve in time. The scope includes algorithm devel- 
opment, assessment of the numerical performance and experimental validation. 
The finite-difference method of discretization is used for the numerical solution of 
differential- equations, while liquid crystal thermography and hot-wire anemome- . 
try are employed in the experiments. Air is considered as the working fluid for 
all the analysis. 

The conjugate gradient method relies upon iterative updating of the unknown 
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functions or parameters through the minimization of a functional. The functional 
is generated in terms of the difference between the computed and measured data. 
Adjoint and sensitivity problems are formulated as intermediate steps to compute 
the gradient of the functional and the step size. The adjoint and sensitivity 
problems also expedite the physical understanding of the problem by enforcing the 
mirror image of the natural law and by estimating the influence of the prescribed 
data on the unknown function or parameter. 

The effectiveness of any inverse technique is considerably dependent upon the 
method of stabilization used. Several tools falling in the category of iterative 
regularization are discussed, devised and employed in the present work. 

The results show that the inverse convection technique is capable of estimating 
the surface information if noise in the additional high quality data is less than 
10%. It is also realized that the additional data must be in the regions of high 
sensitivity. Inverse convection technique has thus been found to be a promising 
tool to bring out useful information from the noisy experimental data. 
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Chapter 1 
Introduction 


While solving a set of partial differential equations, the initial and boundary con- 
ditions must be known. In addition the parameters appearing in the equations as 
well as in the initial and boundary conditions must be specified. Many problems 
of practical interest are governed by partial differential equations. Often they 
suffer from the problem of inadequacy of initial and boundary conditions or the 
knowledge of the required set of parameters. The information may be completely 
unavailable or supplied with substantial scatter. A reasonable solution can how- 
ever be derived if additional good quality data of limited size is supplied. In most 
practical cases, these extra data are also corrupted with some noise. The task of 
an inverse analyst is to come up with an acceptable solution from the body of 
low and high quality data. 

When a physical system is investigated experimentally, the signals generated 
by the measurement devices can hardly be free of any noise. The signal quality 
is quantified by the signal-to-noise ratio. If this ratio is less than 1, the experi- 
mental method needs improvement. If it is greater than 10, filtering techniques 
are adequate to suppress the noise. If this ratio is in between 1 to 10, the inverse 
technique proves effective to subdue the effect of noise. Qualitatively, the scatter 
in raw data can be divided into the following groups: systemic, operational, bias, 
random and sporadic. Poor design of experiment and human factors are respon- 
sible for the systemic and operational scatter, allowing no room to filter it out. 
The bias, random and sporadic scatter are generated from several uncertainties 
associated with the experiment. Fortunately, these types of scatter follow a sta- 
tistically predictive behavior, authorizing some mathematical tools, such as the 
inverse technique, to dampen their corruptive effect. 
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Inverse calculations for thermal problems are of great practical importance 
for scientists and engineers. In thermal sciences, the inverse approach is largely 
applied to. determine thermophysical properties and surface heat fluxes. The 
identification of robust and accurate methods for inverse thermal problems started 
in the 1960s. Different methods emerged, each with built-in advantages and 
disadvantages. In the last decade, attempts have been made to generalize these 
methods. A new and powerful technique based on the conjugate gradient method 
has emerged as a promising tool for inverse thermal problems. The conjugate 
gradient method has several disadvantages aa well. They can be overcome by 
studying the performance of the algorithm in the context of a variety of complex 
applications. 

Inverse thermal problems can be divided broadly into two categories: (1) Pa- 
rameter estimation and (2) Function estimation. In (1), the aim is to determine 
an unknown parameter in the governing equation and/or initial and boundary 
conditions. The determination of initial and boundary conditions are generally 
posed as a function estimation problem. In addition, the determination of ther- 
mophysical properties such as thermal conductivity and heat capacity as functions 
of temperature requires a formulation for function estimation. Occasionally there- 
is an overlap. If a-priori information is available regarding the functional form 
of the boundary conditions, the function estimation problem can be posed as 
a parameter estimation problem. The present work deals with inverse function 
estimation problem and all the discussions and formulations are confined to this 
category. 


1.1 Literature Review 

A survey of the available literature reveals that various aspect of the inverse 
technique have been addressed. These has been presented below as per the fol- 
lowing sections; (1) Mathematical theory of inverse techniques, (2) Inverse heat 
conduction, (3) Inverse convection techniques. 

1.1.1 Mathematical Theory of Inverse Techniques 

The fundamental objective of an inverse procedure is to generate a solution with 
acceptable accuracy and stability from noisy experimental data. Hadamard [1923] 
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introduced the concept of a well-posed problem, originating from the philosophy 
that the mathematical model of a physical problem has to have the properties 
of uniqueness, existence and stability of the solution. Inverse problems generally 
violate these requirements. Often, existence and uniqueness can be forced by 
enlarging or reducing the solution space. For restoring stability, however, one has 
to change the topology of the spaces, which in many cases cannot be enforced due 
to the presence of scatter in the experimental data. The method of regularization, 
introduced by Tikonov [1963] is a watershed in the history of inverse techniques. 
Extensive review of different regularization methods can be found in the work of 
Engl et al. [1990], Bakushinsky et al. [1990], Kirsch [1996] and Glasko [2000] . 

Iterative regularization is the best available option at present as an inverse 
estimation technique. It is based on guessing the function to be estimated and 
improving it by matching the corresponding flow and thermal flelds with exper- 
imental data. In other words the inverse problem can be posed as an optimiza- 
tion, specifically the minimization, of a functional. The very nature of the inverse 
algorithm suggests that small scatter in the measurement can cause enormous, 
perhaps unbounded error in the inverse solution. The method of regularization es- 
sentially adds an additional parameter to the functional to force a stable solution. 
While the regularization method has a firm mathematical background [Kirsch, 
1996, Engl et al., 1990, Bakushinsky et al., 1990] the choice of the regularization 
parameter still needs to be finalized by trial and error. Iterative methods on the 
other hand, always add an error term with the functional and gradually reduce 
the overall error. Iterative regularization stops the iteration at some level where 
this residual error induces stability of the algorithm. However, in this method the 
question of where to stop the iteration, is not clearly understood. The most gen- 
eral approach is to assume equal standard deviation for all the measured data in 
the sense that at convergence, absolute error in each data is less than or equal to 
the standard deviation. While the assumptions regarding the stopping criterion 
of iterative regularization ( also known as discrepancy principle) is not universally 
valid, it is quite popular due to ease of implementation and versatility. 

Steepest descent and conjugate gradient methods are commonly used for iter- 
ative regularization. Different aspects of these methods are extensively discussed 
by Gilyazov et al. [2000]. The conjugate gradient method has the property of 
rapid convergence and a smaller probability of being confined to local minima. 
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This method is now the most accepted for the solution of inverse problems. Hanke 
[1995] has addressed the mathematical basis of the conjugate gradient method, 
the discrepancy principle and their applicability to ill-posed problems. 

1.1.2 Inverse Heat Conduction 

Inverse heat conduction problems are generally formulated to find the surface 
heat flux or thermophysical properties of a solid from interior temperature mea- 
surements. Mathematical models for inverse heat conduction are standardized 
and well-discussed [Alifanov, 1979, Beck et ai, 1985, Kurpisz et ai, 1995, Ozisik 
et al, 2000]. They are also extensively in use. 

Alifanov [1979] was the first to formulate the inverse heat transfer equations 
using the conjugate gradient method. Jarny et al. [1991] outlined the general 
optimization method for solving multidimensional inverse heat conduction equa- 
tions. The idea was to estimate the unknown initial or boundary or conditions 
from internal measurements. They used the conjugate gradient method with the 
adjoint equation. The conjugate gradient method proved to be very efficient even 
without any a-priori assumption about the form of the unknown function. An 
a-priori assumption about the unknown functional serves as a special implemen- 
tation of their formulation. 

Huang et al. [1999] used the conjugate gradient formulation to solve the 
unsteady heat conduction equation in three dimensions. Using the formulation 
of Jarny et al [1991] they successfully estimated the surface heat flux. This w^ 
the first published general solution of inverse heat conduction using the conjugate 
gradient method. 

It is to be realized that the conjugate gradient method is not flawless. One 
major difficulty is to apply perfect stopping criterion of the iterations. Some 
discussions on this topic have been carried out by Khachfe et al [2000]. They 
have pointed out that the minimum value of the functional that will ever be 
attained depends on the amount of error and the initial guess if the data error is 
substantial. 

Khaliday [1998] used a space marching algorithm for inverse heat conduction. 
In the algorithm, stability is achieved by smoothing the measurement data using 
Gram orthogonal polynomials. 
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Lin et al. [2002] successfully coupled inverse technique and transient liquid 
crystal experiment to find the convective heat transfer co-efficient over a sur- 
face. They solved the inverse conduction equation using the Newton root-finding 
technique, ignoring the need for stabilization. 

1.1.3 Inverse Convection Techniques 

Moutsoglou [1989,1990] published the first ever work on inverse convection. She 
applied the function specification method of Beck [1962] in the whole domain as 
well as in a sequential manner. 

Huang et al. [1992] used the conjugate gradient method to determine the 
unknown wall heat fiux in laminar flow with steady forced convection through a 
parallel plate duct. They solved the energy equation for fully developed velocity 
field and developing temperature field. Interior temperature data, very close to 
the plate, were used to construct the functional and to retrieve the boundary 
condition, the unknown wall heat flux in this study. They further simplified the 
problem neglecting the streamwise second derivative term in the energy equation. 
Nevertheless, this simplified approach is the first published work on inverse con- 
vection using the conjugate gradient method. Their study proved the utility of 
the conjugate gradient method in solving inverse convection equations. The au- 
thors also pointed out that conjugate gradient method fails to improve the initial 
guess at some points where the adjoint function becomes zero and proposed a 
modified conjugate gradient approach to circumvent this difficulty. 

Huang et al. [1999] solved a set of unsteady three dimensional inverse forced' 
convection equations and estimated the time varying surface heat flux by conju- 
gate gradient approach. They used commercial software to solve the unsteady 
momentum, energy, adjoint and sensitivity equations in this context. 

Li et al. [2000] applied the inverse forced convection technique to estimate the 
time varying surface heat flux in an annular duct. The energy equation solved 
here was taken to be parabolic in both time and space. The authors used sim- 
ulated temperature measurements at the outer wall of the duct and successfully 
estimated unknown wall heat flux using the conjugate gradient method. • 

Colaco et al. [2001] used a two dimensional inverse forced convection approach 
for estimating boundary heat fluxes in irregularly shaped channels. They used 
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elliptic grid generation procedure for the irregular geometry and determined two 
boundary heat fluxes simultaneously using the conjugate gradient method. The 
solutions were accomplished for time-dependent, spatially dependent and both 
time and spatially dependent heat fluxes. 

1.2 Scope of the Present Work 

In summarizing the previous work on inverse convection technique, it has been 
realized that a few test problems have been addressed and successfully modeled. 
The critical evaluation of the conjugate gradient approach and the inverse con- 
vection technique as a whole, is yet to be presented. In all the previous work 
discrepancy principle is considered as the stabilization method; other techniques 
have remained unaddressed. Furthermore, previous work have used only numeri- 
cally generated data perturbed with random noise as the simulated experimental 
measurements. It is too early to draw conclusions about the performance of in- 
verse convection approach in absorbing the noisy data from real life experiments, 
for which the inverse technique is ultimately called for. Also the potential of the 
sensitivity analysis, in the context of inverse convection is yet to be exploited. 

The present work is concerned with the extraction of useful information from 
raw experimental data using the inverse technique. In this work, the performance 
of the conjugate gradient approach and several stabilization techniques have been 
extensively studied. The inverse algorithm has been applied to simulated as well 
as experimental data. Sensitivity analysis has been carried out to substantiate 
the results. All calculations pertaining to the solution of the partial diflierential 
equations have been numerically carried out. 

With simulated data, the algorithm has been used for estimating the local 
Nusselt number for flow over fiat and ribbed surfaces under conditions of forced 
convection. The idea here is to find the local Nusselt number distribution when 
on one hand the surface temperature carries substantial noise, but additional high 
quality but limited data is available within the flow field. In the context of simu- 
lation, high quality data is generated from the solution of the direct problem with 
an error-free plate temperature. The following test cases have been considered in 
the present study. The direct numerical solvers required at intermediate stages 
of the inverse calculation are also indicated. 
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1. Flow and heat transfer over a flat plate; 

(a) Two dimensional parabolic solver for velocity and temperature; 

(b) Two dimensional flow fleld with three dimensional temperature field 
both determined by a parabolic solver; 

(c) Two dimensional flow field with a parabolic solver with three dimen- 
sional elliptic solver for temperature. 

2. Flow and heat transfer over a flat plate with a single rib: 

(a) Two dimensional elliptic solver for velocity and temperature; 

(b) Two dimensional flow with three dimensional temperature both by 
elliptic solver. . 

(c) Two dimensional elliptic solver for velocity and temperature with tran- 
sient plate temperature. 

In each case the velocity field is considered to be two-dimensional. For fiat 
plate geometry, velocity field is solved through the boundary-layer approximation. 
Flow over a surface mounted rib can only be solved in an elliptic framework. The 
elliptic flow solver is based on stream function- vorticity approach. 

In the experimental part of the study, an inverse technique has been presented 
to estimate the local Nusselt number distribution for flat and ribbed surfaces in 
the laminar flow. The transient temperature field over the surface has been 
imaged by liquid crystal thermography. Two different Reynolds numbers for a 
plate with a solid rib have been considered. 

1.3 Thesis Organization 

The present thesis has been organized in the following manner: 

1. Chapter 1 presents introduction, literature review and scope of the present 
research. 

2. Chapter 2 describes the inverse model in conceptual terms. 

3. Chapter 3 presents the formulation of inverse heat transfer from a flat plate. 
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4. Chapter 4 presents the formulation of inverse heat transfer from a ribbed 
surface. 

5. Chapter 5 is a description of the experimental apparatus. 

6. Chapter 6 presents the experimental data reduction technique. 

7. Chapter 7 reports major results for heat transfer from the flat plate, with 
simulated data. 

8. Chapter 8 reports major results for heat transfer from the ribbed surface. 

9. Chapter 9 carries conclusions of this study and scope of future work. 

10. Appendices contain details of numerical treatment and discretization, code 
validation and grid independence. 



Chapter 2 

Mathematical Formulation: an 
Outline 


The fundamentals of inverse problem and the mathematical tools used to solve 
the resulting equations are described in this chapter. The application of the 
mathematical formulation to specific geometries will be presented in subsequent 
chapters. 

2.1 Fundamentals of Inverse Problem 

Following Keller [1976], two problems are said to be inverse to each other if the 
formulation of one involves the other. For mostly historic reasons, one might 
call one of these problems (usually the simpler one or the one that was studied 
earlier) the direct problem; the other is the inverse problem. However, in the real- 
world, there lies, in most cases, a natural distinction between the direct and the 
inverse problems. If one wants to predict the future behavior of a physical system 
from knowledge of its present state and physical laws (including concrete values 
of all relevant parameters), it will be called a direct problem. Possible inverse 
problems are the determination of the past state from present-day observations 
(i. e., the calculation of the evolution of the system backward in time) or the 
identification of physical parameters from observations of the evolution of the 
system (parameter identification). Thus inverse problems are associated with 
determining the cause from the desired or an observed effect. 

The detailed classification of inverse problems is given by Kozdoba[1982]. In 
general, they fall into the following categories: 
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1. Inverse Boundary Value Problems. 

2. Parameter Estimation. 

3. Inverse Initial Value Problems. 

4. Inverse Geometry Problems. 

5. Others. 

In the present work the inverse problem is formulated as an inverse boundary 
value problem. Inverse boundary value problems can be viewed as the recovery of 
boundary conditions of a partial or ordinary differential equation from complete 
or partial interior information. In heat transfer, this generally refers to estimation 
of surface heat flux history of a body from transient temperature measurements 
at some interior locations. 

2.2 Ill-posedness of inverse problem 

Hadamard [1923] introduced the concept of a well-posed problem, originating 
from the philosophy that, the mathematical model of a physical reality has to 
have the properties of uniqueness, existence and stability of the solution. If one 
of the properties fails to hold, the problem was called ill-posed. Let us define a 
problem by the operator equation 

Db = e (2.1) 

where b E B and e E E, E and B are metric spaces and D is an operator 
so that DB = E. In general, b can be a vector that characterizes a model of a 
phenomenon and e can be the observed attributes of the phenomenon. A correctly 
posed problem must meet the following requirements: 

1. The solution of Equation 2.1 must exist for any e E E. 

2. The solution of Equation 2.1 must be unique. 

3. The solution of Equation 2.1 must be stable with respect to perturbations 
on the right-hand side, i.e. the operator must be deflned in the space 
E and be continuous. 
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Let us assume that the values of e cannot be exactly measured and only the 
approximate values of 6 E. If the error is bounded 

\\e-es\\E<5 (2.2) 

where ||...|| is the Euclidian norm, defined in a metric space E. If the problem is 
correctly posed the following condition holds: 

I16-5j||£-^ 0 if 5 ^0 (2.3) 

where b is the solution of equation Equation 2.1 for an accurate vector e and 

bs = D-^s (2.4) 

This result (Equation 2.3) is not true when ill-posed problems are encountered. 
For ill-posed problems the inverse operator D~^ defined throughout its domain 
DB e E is not continuous, it means that the solution of Equation 2.1 does not 
depend continuously on the data. Thus, the element bs = D~^es, even if it is 
exists, is not an approximate solution, because bs may arbitrarily deviate from 
the exact solution b even for small S. 

A lot of problems of optimum control, linear algebra, integral equations, the 
summation of Fourier series with approximately specified coefficients and the min- 
imization of functionals are ill-posed and they are extremely sensitive to para- 
metric and measurement errors. 

The ill-posed nature in a mathematical sense of inverse problems originate 
from physics. Two effects should be mentioned here: damping and lagging. In 
the direct problems of heat transfer, the temperature fiuctuations are diminished 
internally in comparison with the surface temperature changes. This is a damping 
effect. In particular, the higher frequency components of the boundary conditions 
are damped at a higher rate than the lower frequency components. In the inverse 
problem the higher frequency components of the measurements are amplified as 
the measurements are projected to the surface. As a result, the surface condition 
predictions might be overwhelmed by the noise in the measurements. 

A large time delay or lag in the internal response to the changes in the bound- 
ary conditions should a.lso be noticed. As an example, let us consider the problem 
of heating a symmetrical slab of thickness 21. The Biot number is defined as 
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^1 = hi/ k and the Fourier number is defined as Fo = at / where h is the heat 
transfer coefficient, t represents time, a stands for the thermal diffusivity and k is 
thermal conductivity of slab material. It can be proved that if the boundary con- 
dition of the third kind is prescribed on the surface of the slab, the temperature 
on the axis of the slab remains almost unchanged until Fo Ri 0.05. At the same 
time the surface temperature for example, Bi=4 becomes equal to one half of the 
temperature for the steady-state. As a result the surface temperature cannot be 
estimated from the temperature measurements on the axis of the body within 
such a time interval. 

Another question that arises in inverse thermal problems concerns the unique- 
ness of surface condition as predicted by discrete internal measurements. The 
natural way to overcome this difficulty is to obtain the maximum amount of 
information from the measurements. For transient problems, this implies mini- 
mizing time steps between two subsequent measurements. However, small steps 
may cause instabilities in the two subsequent measurements, causing instabilities 
in the solution. This observations is contrary to the direct problem where sta- 
bility can be improved by reducing the time steps. As a consequence, for inverse 
problems the upper and the lower limits for time steps need to be specified. 

The ill-posed nature of an inverse formulation causes elegant methods of di- 
rect problems be inapplicable. Special numerical techniques must be employed 
to stabilize the calculations. Any procedure for inverse problems, given in the 
literature without discussion on the question of stability is, in fact, without value. 

2.3 Methods of Solution 

Since numerous inverse methods have been proposed in the literature, some cri- 
teria for their evaluation are needed. The criteria listed are those proposed by 
Beck [1985]: 

1. The predicted temperatures and heat fluxes should be accurate if measured 
data is of high accuracy. 

2. The method should be insensitive to measurement errors. 

3. The method should be stable for small time steps or intervals. This permits 
the extraction of more information regarding the time variation of surface 
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conditions than is permitted by large time steps. 

4. Temperature measurements from one or more sensors should be permitted. 

5. The method should not require continuous first order time derivatives of the 
surface heat flux. Furthermore, step changes or even more abrupt changes 
in the surface heat fluxes should be permitted. 

6. Knowledge of the precise starting time of the application of the surface heat 
flux should not be required. 

7. The method should not be restricted to any fixed number of observations. 

8. Composite solids should be permitted. 

9. Temperature dependent properties should be permitted. 

10. Contact conductances should not be excluded. 

11. The method should facilitate easy computer programming. 

12. The computer cost should be moderate. 

13. The user should not have to be highly skilled in mathematics in order to 
use the method. 

14. The method should be capable of treating various multidimensional coor- 
dinate systems. 

15. The method should permit extension to more than one heating surfaces. 

16. The method should have a statistical basis and permit various statistical 
assumptions for the measurement errors. 

Though the above criteria are proposed for inverse conduction problems, ma- 
jority of these are applicable for inverse convection as well. The formulation used 
in the present work meets most of these conditions (with the exception of (12)). 

While solving direct thermal problem one has to formulate and analyze gov- 
erning partial differential equations of the related phenomenon. These equations 
must be supplemented with appropriate initial and boundary conditions. Ail 
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coefficients in the governing equations, and initial and boundary conditions in- 
cluding thermal conductivity, density and specific heat as well as the capacity 
and the location of internal heat sources, if they exist, have to be specified. Also, 
geometry of the physical domain must be fully specified. The main objective 
here is to find the temperature distribution within a domain initial and boundary 
disturbances. 

For inverse problems, some of the above listed input data are missing. Thus, 
the mathematical description of the direct problem becomes incomplete and defies 
its solution. To make the solution of the inverse problem available, additional 
information is required. Usually, this data is obtained from measurements at 
interior points of a domain. The general strategy of solving inverse problems can 
be summarized by the following steps: 

1. Make the mathematical description of the boundary value problem complete 
by assuming arbitrary values as required by the direct problem but not 
specified in the input data of the inverse problem. 

2. Solve the direct problem by an analytical or a numerical method. 

3. Compare the calculated and measured values and modify the assumed input 
data to ensure the best match between these quantities. This essentially 
calls for an optimization scheme to minimize the functional created from 
the difference of the computed and the measured data. 

4. Stop iteration at a level where maximum information and minimum noise 
is passed to the numerical solution. 

Following this strategy of solving inverse thermal problems, one not only dis- 
covers the cause from a known result but also retrieves missing information. 
Particular methods of solution of inverse problems differ way the assumed input 
data is mo dified^TFhe techniques for-so lving^n inverse problem can^he loosely 
classified under the following groups: 

1. Integral equation approach. 

2. Integral transformation techniques. 

3. Series solution approach. 
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4. Polynomial approach. 

5. Hyperbolization of the heat conduction equation. 

6. Space marching technique with filtering of the noisy input data, such as the 
mollification method. 

7. Iterative filtering techniques. 

8. Steady-state techniques. 

9. Beck’s sequential function specification method. 

10. Levenberg-Marquardt method for the minimization of the least-square norm. 

11. Tikonov’s regularization approach. 

12. Iterative regularization methods for parameter and function estimations. 

13. Neural networks and genetic algorithms. 

The details review of these methods can be found in the works of Alifanov 
[1979], Beck et al. [1985], Kurpisz [1995] and Ozisik et al. [2000]. Presently, the 
iterative regularization method has become popular among researchers, due to 
its simplicity and versatility. In the present work iterative regularization method 
has been used with some traditional and modified schemes of stabilization. 

2.3.1 Iterative Regularization 

The iterative regularization method, developed by Alifanov [1979], essentially 
converts the inverse problem into a statement of optimization problem. It is solved 
through a series of direct problems and a suitable stopping criteria. For well posed 
problems the iterations are carried out as long as the discrepancy between two 
successive steps becomes sufficiently small. For ill-posed problems the number 
of iterations P must be coordinated with the data error <5. In that sense, the 
number P acts as a regularization parameter. The iterative regularization can be 
applied to inverse boundary value problems as well as others. 

Let the problem to be solved be expressed by an operator equation 


Dh = e^ 


(2.5) 
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If the operator D ^ is not continuous then the problem becomes ill-posed. To 
overcome the ill-posed nature of the problem two questions have to be answered: 

1. What iterative algorithm should be employed? 

2. How many iteration steps should be performed? 

Numerous iterative algorithms [Kurpisz, 1995] can be applied to solve Equa- 
tion 2.5 for b. Some of them are listed below: 

1 ^+1 _ (DbP -es),b<^ = 65 , p= 1,2,... P 

where P stands for the total number of iterations and p refers to the current 
iteration level. 

2. = (I - uD*D)lf + vD*es, b° = D*e 5 ,p =1,2,...P 

where D stands here for a matrix operator, 6 is a vector, I is an identity 
matrix, D* is the conjugate of matrix D, v is a. coefficient that satisfies 

0<u< 2/l|T>*D|| 

3. If*' = 5ii ELoK^ - + ue,] 

The problem of solving an operator equation can be transformed into the 
problem of minimizing the functional 

J=\\Db5-e5\\ (2.6) 

It leads to the successful iterative methods, called the gradient search techniques. 
The general form of an iterative algorithm can be written as follows: 

6?+^ = 5) (2.7) 

The iterative algorithm is termed as regularized if it satisfies the following 
conditions: 

1. For any 5° G B the iterative sequence tends to the exact solution when- 
S = 0 

2. The transition operator i?(5?, 5) is continuous for each element b when 5 = 0 
These conditions are satisfied by the following well-known methods: 
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1. Steepest Descent. 

2. Conjugate Gradient. 

Application of the steepest descent method results in the following iterative 
formula: 


= If - PfJl, p = 0,l.. 

.P 

(2.8) 

e _ II-CIII 

’ nojni 


(2.9) 

The conjugate gradient method produces the following algorithm: 


tf*^ = if - p = o,i.. 

.p 

(2.10) 

21|DQP||| 


(2.11) 

= ao = 0, = 

IW'III, 

(2.12) 

ii.'r'iii 


where 

J',^ = 2D*{D}fs-es) 

is a gradient of the functional 

jr=\m-es\\ 

with respect to b if the operator D is linear. The symbol D* stands for the 
conjugate operator, whereas the notation (a, b) stands for a scalar product. Sim- 
ilar formulas can be written for cases when the operator D is nonlinear but the 
convergence of the iterations must be established via numerical experiments. 

As indicated, the number of iteration steps P must correspond with input 
data errors and to evaluate its value the discrepancy principle can be used. This 
means that to determine F we have to take the first solution at the P the iteration 
level that satisfies 

- e^ll < 0:5 (2.13) 

where a > I is a fixed constant. 
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2.4 Conjugate Gradient Approach 

The gradient of a function is a very important property. The partial derivative of 
a function / with respect to each of the n independent variables are collectively 
called the gradient of a function and is denoted by v/- If we move along the gra- 
dient direction from any point in n-dimensional space, the function value changes 
at the greatest rate. Hence the direction is called the direction of steepest ascent. 
Since the gradient vector represents the direction of steepest ascent, the negative 
of gradient vector denotes the direction of steepest descent. Thus any method 
that makes use of the gradient vector as a descent direction can be expected to 
yield rapid convergence. The steepest descent method for minimization of a func- 
tional is based on this criterion. Unfortunately, the direction of steepest descent is 
a local property, and not a global one. Many improvements have been suggested 
to the steepest descent method. One of them is to consider both the present and 
the previous descent directions to find the new descent direction. This technique 
is known as the conjugate gradient method. In the present work this technique 
is used to solve the inverse convection equation. The conjugate gradient method 
has been used to optimize the functional formed by the difference between mea- 
sured and computed temperatures and the governing partial differential equation 
with the boundary and/or initial conditions serving as the constraints. The entire 
algorithm is divided into a number of subproblems as described below. 

2.4.1 Direct Problem 

The direct problem is to solve the governing partial differential equation with the 
proper boundary and/or initial conditions. In the context of an Inverse calcula- 
tion, guessing the missing informations, which are continuously updated. 

2.4.2 Inverse Problem 

If 9 and Y are the computed and measured temperatures respectively the follow- 
ing functional can be formed: 

[{e-Y)H{v-yi) (2.14) 

Jxi 

where, 5 is the Dirac delta function and yi indicates the plane of measurement. 
The integration is carried out over this plane. The solution is generated by seeking 
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the function 9{x, y) that minimizing J. 

2.4.3 Sensitivity Problem 

The sensitivity problem is obtained by perturbing the governing equation and the 
initial and boundary conditions of the direct problem. The sensitivity function A0 
is the solution of sensitivity problem, and is defined as the directional derivatives 
of the temperature 6 in the direction of perturbation. In the algorithm, the 
sensitivity function is used only to determine the step size. The solution procedure 
of the sensitivity problem is similar to that of the direct problem, but the physical 
implication of the sensitivity function is enormous. This function dictates the 
dependence of temperature on the unknown quantity to be estimated (here, the 
local Nusselt number). Specially in thermally developing flows, the importance of 
the sensitivity function is revealed because temperature within only a very small 
region of the domain is influenced by the local Nusselt number distribution. The 
sensitivity function discovers the interior region where the boundary information 
is contained. In this way, it helps in designing the experiment for high quality 
measurement. The sensitivity function also indicates the influence of data error 
on the final solution and helps in selecting the best method of stabilization. 

2.4.4 Adjoint Problem and Gradient Equation 

The formulation of the adjoint problem proceeds by using the principle of vari- 
ational calculus. Here the concept of Lagrange multiplier is introduced. The 
Lagrange multiplier is obtained by solving the the adjoint equation with initial 
and boundary conditions. The Lagrange multiplier serves as the input to the gra- 
dient equation to evaluate the gradient of the functional. The gradient in turn is 
used to find the descent direction. 

2.4.5 Step size computation 

The gradient of the functional serves as an input to the sensitivity problem. The 
step size is computed in relation to the sensitivity coefiicient and the differences 
of the measured and computed temperature, with the requirement that the step 
size applied to the gradient would bring the functional to its minimum. 
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2.5 Algorithm Stabilization 

Stabilization is the method of imposing stability to the inverse algorithm, specif- 
ically the flow of scatter into the flnal solution. The aim is to resist unbounded 
changes in the solution for a small change in the input data. Clearly, if the exper- 
imental data is completely free of error, no stabilization procedure is required and 
minimization of the functional would produce the correct result. Unfortunately, 
in practice, experimental data can hardly be exact. The idea of regularization 
involves (1) modifying the functional before carrying out optimization or (2) pre- 
mature stopping of the minimization process. A few approaches used in the 
present work are discussed below. 

2.5.1 Discrepancy Principle 

The discrepancy principle is used by both Tikonov [1963] and Alifanov [1979]. 
Presently this is the most widely used method of stabilization. The principle is 
based on the following assumptions; 

1. Standard deviation (cr) of all measurements are equal. 

2, The value of this standard deviation is known. 

The iterative procedure is stopped when the following criterion is satisfled 

(2.15) 

The limitations of this procedure comes from the assumptions. Data scatter * 
in different portion of the physical region may not be less than the standard 
deviation. Even if they are comparable, their standard deviation may not be 
known in advance. 

2.5.2 Additional Measurement Approach 

This method proposed by Ozisik [2000] circumvents the limitations of the dis- 
crepancy principle. Here two sets of data are used for constructing functionals Ji 
and J 2 . While minimization of Ji is continued, the value of J 2 initially decreases. 
At some point the value of J 2 starts to increase, revealing the onset of instability. 
Iteration is generally stooped at this stage. ' 
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While this method has a number of advantages over the traditional discrep- 
ancy principle, it needs a large number of data. It can fail to give satisfactory 
results if the orders of magnitude of data scatter are different. 

2.5.3 Smoothing of Measurement Data 

The stability of inverse algorithm is hampered by the high frequency components 
of the measured data. If the high frequency components can be cut out to some 
extent, the solution can be stable. This idea has been exploited by Khaliday 
[1998] for inverse conduction. 

2.5.4 Present Approach 

In the present work, two types of data has been handled: numerically simulated 
data and experimentally generated data. Different methodologies are employed 
for these data types. 

Numerically Generated Data 

Here the standard deviation of the data set are not all the same. A large 
number of data with high error, whereas number of high quality data is less. If 
we use the stopping criteria described above the solution would be overwhelmed 
by the data with high error and small number of high quality data would fail to 
refresh the solution. To circumvent this difficulty, the following definition is used 
to formulate the functional. 

J = / [m(« - YfSiv - m) +{6- Y)H(v - !/2)) (2.16) 

Jn 

where w is a weight assigned to the high quality data. 

The data at y = y 2 is taken to contain substantial error while the data at 
y = yi is of high quality. Before constructing the functional all data is filtered 
using Gram orthogonal polynomial filtering. With due weightage through the pa- 
rameter w and filtering, the functional is minimized to yield the inverse solution. 

The determination of the parameter w is very crucial. This determines the 
effect of different levels of error on the final solution. The method of minimization 
of the functional through the conjugate gradient method has shown that the 
minimum value of the functional {Jmin) depends upon the data scatter (Figure 



2.5 Algorithm Stabilization 


22 


2.1). This phenomenon was addressed by Khachfe et al. [2000]. Jmin reaches a 
value close to zero if the scatter is negligibly small. As the data scatter increases 
the absolute value of Jmin increases. Though the mathematical basis of this 
observation is yet to be addressed, this property is exploited in this present work 
to evolve an idea of the value of w. 



Figure 2.1: The variation of functional with iterations as shown by Khachfe et 
al. [2000] 


Starting from an arbitrary initial guess, let the minimum value of the func- 
tional be Jimin, when data with high scatter is used and J 2 mtn) When only data 
with low scatter is used. Allowing both the data set to influence the flnal solution 
equally, we can use 

'^optimum ^ "T" (2*17) 


J Imin 

If w « m optimum the noisy data will influence the solution very much and 
there will be little recovery of the local Nusselt number. If w »• w optimum 
the noisy data is practically kept out of any use and the solution is based on 
very few high quality data points. This set might not be capable of producing 
uniqueness to the solution. Observations of the present work show that the range 
.y; = (10 - 100)woptimum gives the best estimations with accuracy and uniqueness 
both being reasonably maintained. 
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Experimentally Generated Data 

For experimentally generated data the standard deviation of the data is not 
clearly known. A large amount of data is available. Additional Measurement 
Approach, described in Section 2.5.2, has been used here. 

2.6 Filtering of Data 

The objective of filtering (or, smoothing) is to cut out the high frequency compo- 
nents from the measured data. The high frequency components accentuate the 
errors in the solution of the inverse problem to a great extent. As shown in the 
work of N. A. Khalidy [1998] and also of J. Taler et al. [1999], the filtered data 
does stabilize the algorithm. The present work also supports this observation. 

Several filtering techniques are presented in the literature to eliminate the 
noise from the measured data. Examples are inclusion of digital filters in the 
inverse algorithm. Hanning, triangular, Gaussian, Tukey, Savitzky-Golay and 
Gram orthogonal polynomial are examples of such filters. In the present work, 
the performance of Hanning and Gram orthogonal polynomial filtering have been 
tested and the later is chosen. 

2.6.1 Hanning Filter 

Hanning filter uses the concept of frequency domain filtering using fast Fourier 
transform. This is extremely suitable for cases where the number of data points is 
large. A stringent requirement of this filtering technique is that, for using the fast 
Fourier transform algorithm, it needs number of data to be equal to 2”, where n 
is a positive integer. The fast Fourier transform transfers the data into frequency 
domain. After cutting out the highest frequency portion (approximately 10% of 
the area under the power spectrum), the data is again transfered to time or space 
domain by means of inverse Fourier transform to get the smoothed data. This 
method has other difficulties such as aliasing and distortion of phase information, 
and should be carefully employed. 
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2.6.2 Gram Orthogonal Polynomial Filter 


Least squares approxiraation is very well suited for the recovery of a smooth 
function from noisy information. It is possible to chose an appropriate function 
which is flexible enough to reconstruct the underlying noise-free function and 
its derivatives while still being orthogonal to the noise, i.e., fails to follow the 
oscillations in the measured data. The Gram orthogonal polynomial filtering, a 
popular methodology of stabilization of inverse problem [Khalidy, 1998], follows 
the principle. It uses weighted averaging of several past and future (spatially or 
temporally) data for smoothing. Thus the filter replaces each data value Tj by a 
combination of itself and certain adjacent points by fitting a piecewise polynomial 
using the least squares method to find the coefficients of the polynomial. It can 
be shown that, 

n 

= ( 2 . 18 ) 

i=o 


where, fnixi) is the smooth representation of the measured temperature at a; = Xj. 
Let L be the number of points used to the left and to the right of the central 
point X = Xi. (Notice that Equation 2.16 can be used for smoothing the data 
spatially as well as in time.) The subscript n refers to the total number of data 
points used for the smoothing, the process (n = 5 for L = 2, n = 7 for L = 3, and 
so on). The parameters in Equation 2.16 can be calculated as [Korn, 1968], 


b' - ( 2 Z,+ 3 +l)!( 2 n- 3 )! 
( 23 + 1 ){( 2 L )!}2 


and. 




i)[2’”^(L -h 2)N 
(m!)2(2L)N 


(2.19) 


.( 2 . 20 ) 


where = y(y — l)(y ~ 2). {y — m -f- 1) and m — 1, 2 This type 

of smoothing is also called moving window averaging. Choice of number of data 
points is generally heuristic; 7 to 11 points are expected to yield reasonably good 
approximation. 

Gram orthogonal filtering by 7 points moving window averaging is employed 
in the present study. In this method a piecewise polynomial is set through the 
first 7 points and these 7 points are smoothed. The next 7 points are then taken 
into account. It is assumed that at the beginning and the end the data points 
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are fixed. Using the above Equations 2.16, 2.17 and 2.18 the smoothed values of 
fn can be computed as: 

fnixi-s) = ^(39Ti_3 + 8Ti_2 - 4ri_i - 4Ti + + iTi+2 - 2Ti+,) (2.21) 

fnixi-2) = ^(8Ti_3 + 19Ti_2 + 16Ti_i + QTi - 4Ti+i - 7Ti+2 + (2.22) 

fnixi.i) = ^(-4Ti_3 + 16Ti_2 + 19Ti_i + 12ri + 2Ti+i - 4Ti+2 + T^+a) (2.23) 
fnixi) = — (— 2Ti_3 + 3Tj_2 + 6Ti_i +• 7Tj + QTi^i + 3Tj+2 — 2Ti^3) (2.24) 

fn(xi+i) — —{Ti-3 — 4rj_2 + 2Ti_i + 12Ti + 19Tj^.i + 16Tj+2 “ 4T’i^.3) (2.25) 

fn{xi+2) = “ 7Ti_2 — ^Ti-i + 63i + 16Ti+i + 19Ti+2 + (2.26) 

fnixi+s) = — (— 2Tj_3 + 4ri_2 + Ti_i — ATi — 4Ti4.i + 8Tj+2 + 39ri+3) (2.27) 

Near the end of the domain, the number of points available may not be suf- 
ficient to carry out this formulation. Reasonable results can still be obtained by 
extrapolating the measured data set[1998] 

2.7 Computational scheme and Discretization 

For the numerical solution of direct, sensitivity and adjoint problem the governing 
partial differential equations and boundary conditions have been discretized using 
finite difference technique. For two and three dimensional parabolic equations the 
Crank-Nicholson scheme is used; fully implicit scheme has been used for elliptic 
problems. For the flat plate, the velocity field is computed using the boundary- 
layer approximation. For the ribbed plate, a vorticity-stream function approach is 
employed for the velocity field using an explicit scheme for the vorticity equation. 
In all cases the flow field is assumed to be laminar and buoyancy forces have been- 
neglected. Steady and unsteady forced convection heat transfer conditions have 
been considered. The numerical solution of the direct problem is discussed in 
Appendix A. 



Chapter 3 


Flow and Heat Transfer over a 
Flat Plate 


The mathematical procedure of the inverse technique based on iterative regular- 
ization method is presented in this chapter for flow over flat plate. The origi- 
nal algorithm, using conjugate-gradient gradient method for conduction problem 
was developed by Alifanov [1979], which was later applied to convection problem 
by Huang et al. [1992, 1999]. The following assumptions have been employed 
throughout: 

1. The fluid is Newtonian and incompressible. 

2. The fluid properties are constant with respect to temperature. 

3. The flow is laminar. 

4. The velocity field is steady and two-dimensional. 

5. The fluid inlet temperature is spatially uniform and temporally constant. 

6. Viscous dissipation and mechanical work terms are negligibly small. 

7. Forced convection is predominant; hence the momentum and energy equa- 
tions are not coupled. 

8. The flat plate temperature is temporally a constant. 

9. Plate temperature data is available with some scatter whose characteris- 
tics are known; in addition, limited but high quality temperature data is 
available within the flow field. 
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The algorithms described in this chapter, consists of the following basic steps: 

1. Smooth the given temperature data using Gram orthogonal polynomial 
filter through the seven point moving window averaging scheme. 

2. Solve the continuity and momentum equations to find the velocity field. 

3. Take an initial guess of local Nusselt number distribution over the plate. 

4. Compute the temperature field using a guessed local Nusselt number as the 
boundary condition at the plate. 

5. Compute the optimization functional. 

6. If the value of the functional conforms to the stopping criteria, stop itera- 
tions and present the results. 

7. Otherwise, solve the adjoint problem and find the gradient of the functional. 

8. Find the search direction. 

9. Solve the sensitivity problem and find the search step size. 

10. Update the local Nusselt number values over the plate surface. 

11. Return to 4 and continue. 


3.1 Velocity Computation 


The boundary layer approximations are assumed to be applicable for the veloc- 
ity field, i.e. the streamwise second derivatives in the momentum equation are 
neglected. This step converts the original elliptic momentum equation into a 
parabolic form, thus enabling a marching solution. However, this restricts the 
investigation to high Reynolds number. The boundary-layer form of continuity 
and momentum equations in the non-dimensional form can be written as: 


du dv 
Bx By 


(3.1) 


B(u^) B{uv) _ 1 d'^u 

Bx By ~ Re' 9 y 2 


(3.2) 
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Figure 3.1: Flow over a Flat Plate 


The boundary conditions with reference to Figure 3.1 are: 


a; = 0, u 

= 1, v = 0 

(3.3) 

y = 0, u - 

= 0, u = 0 

(3.4) 

y = d, 

, ti = 1 

(3.5) 


3.2 Steady Two-Dimensional Temperature Field 

The inverse formulation for a steady thermal field is presented first. 

3.2.1 Assumptions 

1. Temperature field is two-dimensional, i.e, the plate temperature is varies 
only with length x but not z or time. 

2. The temperature field is everywhere steady. 

3. Boundary layer approximations are applicable to the thermal energy equa- 
tion. 
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Figure 3.2: Heat Transfer over a Flat Plate 


3.2.2 Direct Problem 


The direct problem solves the following energy equation: 

d{u6) d{v9) _ 1 d'^d 

dx dy Pe dy"^ 

The boundary conditions with reference to Figure 3.2 are: 


(3.6) 


X = 0, 6 = 0 


(3.7)' 


y = 0, ^ = -Nu (3.8) 

dy 

y = d, 6 = 0 (3.9) 

where Nu is the local Nusselt number. In an overall inverse calculation, local 
Nusselt number values are updated with the iterations. 


3.2.3 Inverse Problem 

The inverse problem is posed in the form of minimization of the following func- 
tional: j 

j = / +w(e„-Y„f]ix , (3.10) 

Jo 

where and Tu, are the computed and measured temperatures at the wall, 6m 
and Tm are the computed and measured temperatures in the interior and to is 
the weight assigned to the high quality data. 
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3.2.4 Sensitivity Problem 

The sensitivity problem determines the the dependence of the temperature on 
the unknown local Nusselt numbers. To formulate the sensitivity problem, a 
perturbation of is given to Q and Q is replaced by ^ -f- in Equation 3.6 to 
yield 

d[ui9 + A9)] d[v{9 + A9)]_ld^{9 + A9) 

dx ^ dy ~ Pe dy^ ^ 

Now, subtracting Equation 3.6 from Equation 3.11, we have the following sensi- 
tivity problem 

a(«A^) ^ 1 a^Ae) 

dx dy Pe dy"^ 

In the boundary condition, Nu is replaced by Nu-t-ANu and 9 by 9 + Ad; thus 
we have: 

x = 0,9 + A9 = 0 (3.13) 

d{9 + A9) _ 


2 / = 0 , 


dy 

y — d, 9 + A9 = b 


(Nu + ANu) 


(3.14) 

(3.15) 


Now, subtracting Equations 3.7, 3.8, 3.9 from Equations 3.13, 3.14, 3.15 respec- 
tively, we get the following boundary conditions for the sensitivity problem: 


a: = 0, A0 = 0 


2 / = 0 , 


y = d, A9 = 0 
d{A9) 


dy 


ANu 


(3.16) 

(3.17) 

(3.18) 


3.2.5 Adjoint Problem 

Let, A be the Lagrange multiplier that is associated with the constraint. We can 
write the functional 

J = r [‘'[i9^ -Y^)^ + w{9m-Ym)^]dxdy + 

Jo Jo 


1 d'^9 d{u9) d{v9) 

A 


If 


[Fedy'^ dx dy J 


dxdy 


(3.19) 
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In Equation 3.19, setting 6 + A6 in place of 9 and J + AJ in place of J, we 


J + A J = r r [{9 + AO- Yf5{y) +w{9 + A9- Y)^5{y - yi)] dxdy 
Jo Jo 

1“^ ^ p 3^9 + A9) _ d{ui9 + Ae)} _ + 

Jo Jo [Pe dy2 dx dy \ 

(3.20) 


where, S{y — yo) is the Dirac delta function centered around yo 

Subtracting Equation 3.19 from Equation 3.20 and neglecting the terms in- 
volving squares of A^, we have 


AJ = 2 


A9 [w{9 - Y)5{y) + {9 - Y)5{y - yi)]dydx + 


Jo Jo [Pe dx dy j 


Using the boundary conditions from sensitivity problem, we get 


Jo L % Jj,=0 
Jo L oy 

= rh^r\. 

Jo L oy Jj,^o 

n d a2\ 

" 5 ?** 




diA9) dX 


dxdy 


lo Jo dy dy-'^' 
f^l^^dX fd'^X 

L . dy J dy^‘ 

r\^A9\ dx + 

lo [dy Jj,=o 


- y=d 

A9dy dx 

- y,==0 


(3.22) 


Further 


^d{uA9) 
) Jo dy 


X — dxdy = / / AuA0 dy 

dy Jo L J J®=o 

1 . r^ dx 


[AuA^]“ J dy - —uA9dxdy 

f.d dX 

[AwAl9]^_i dy - / / —uA9dxdy (3.23) 
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Also 




X=1 


dy 


x=0 


1 f*d 



dX 


10 Jo 9y 

Now combining Equations 3.20 to 3.24, we have 


vAddxdy 


(3.24) 


AJ 


rdr I q2x dx dX 
+ U-;r- + V- 


If 

Jo Jo L 


Pedy"^ dx dyj 


AOdxdy -t- 


2 / [ [w{6 -Y)5{y) + {9 -Y)5{y - yi)]A6dxdy + 
Jo Jo 

Ae^ [- 

Jo Pe [< 


'dX' 


r a(A«)i 


Pc h 



dx + 


[ [>^uAd]^^^dy 
Jo 

p 1 r 1 d'^X d{uX) I d{vX) 

Jo Jo 


+ 


dx 


dy 


+ 2{e-Y)Siy-y,) 


AOdxdy -f- 


/' 


X 


Pe dy 


4 - 2w(9yj — Yu)) 


d{A9) 

dy 

AOdx 


y^d 


dx 


y=0 


(3.25) 


y=0 


where, 9yj and Yyj are computed and measured temperatures at the plate, respec- 
tively. For AJ -4 0, the following adjoint equation and boundary conditions can 
be deduced. 


i.A+®(!N + ^ + 2(«-y)%-yO = 0 

PeSy2 ^ dx dy ^ ^ ’ 

a; = 1, A = 0 

o Y 

y = 0, — = — 2ti;Pe(0ju — Yyj) 
y = d, X = 0 


Since 


AJ= -^ J^ [A(-Aiyu)]^^o <^3: 

and the definition of gradient 

AJ= J'i-Am)dx 

Jo 


(3.26) 

(3.27) 

(3.28) 

(3.29) 

(3.30) 

(3.31) 
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From Equations 3.30 and 3.31, we have the following expression for the gra- 
dient of the functional, 

AJ A(a;,0) 


J' = lim 


(3.32) 


aNu-^o ANu Pe 
The adjoint equation deduced above has negative diffusion. It has no initial 
condition either. These are circumvented by the transformation x — 1 — x, which 
coverts the equation to a well behaved form, generates an initial condition and 
makes the diffusion term, to be positive. 

3.2.6 Determination of Search Step Size. 

If the value of the functional is J* after k-th. iteration and it is after tem- 
peratures are updated through a step size of the condition for the optimum, 
step size will be given by: 


5/5^ 


0 


(3.33) 


Recalling the expression of functional, 

pd 


Jfc = f [' [vj{e - Y)H{y - Vi) + (^ - Y)Siy)] dxdy (3.34) 
Jo Jo 

minpJ>^{-ANuf = minp [' [\w[e{i-ANu)‘^ - -YfSiy - y^) 

Jo Jo 

+[e{{-ANu)‘^ - p’^P'^} - Yf5{y)}dxdy (3.35) 


Linearizing by the Taylor’s series expansion, 

pd 


jfc+i[(_AiVu)*’+^] = w [' [\9[i-ANuf~P'‘A9{P'^)]-Y}H{y-yi)dxdy 

Jo Jo 

+ r / {9[i-ANu)>= -p'=A9{P^)]-Y}^6{y)dxdy 

Jo Jo 

- [ [ {{9 ~Y) ~ p’°A9}^{w6(y - yi) + 6(y)}'dxdy 

Jo Jo 


(3.36) 


2 [ f A9{9 - Y p’‘A9){w5{y - yi) + S(y)}dxdy (3.37) 
Jo Jo 

Using Equation 3.33 and 3.37, we have the following equation for the step 




Size: 


p' 


>k _ - - Vi) + 

J^{A9)^{vj6(y - yi) + 6{y)}dxdy 


(3.38) 
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3.2.7 Complete Inverse Algorithm, 

The complete inverse algorithm for calculation the Nusselt number distribution 
over fiat surface is summarized below: 


1. Solve the continuity and momentum equations to get the velocity field. 

2. Filter the^ imperfect plate temperature data using Gram orthogonal poly- 
nomial filter. 

3. Solve the energy equation using the following boundary conditions: 

a; = 0, 0 = 0 (3.39) 

y = 0, 9 = Oyj (filtered) (3.40) 

y = d, 9 = 0 (3.41) 


4. Get the values of local Nusselt numbers. Take an initial guess of Nusselt 
number based on the distribution obtained. Initialize J = 0 

5. Solve energy equation using the following boundary conditions (with Nusselt 
number from step 4); 

x = 0, 9 = 0 ■ (3.42) 




(3.43) 

y = d, 9 = 0 (3.44) 

Find the temperatures at the plate as well as the plane of high quality data. 

6. Compute the functional 

J = [ [ {{^W - Ywf + iu{9ra - Ymf}dxdy (3.45) 

Jo Jo 


7. Check the stopping criterion 


\r-r-^\<e 

where, e is a predetermined real number, very close to zero. 
If it is satisfied, stop iterations. Otherwise 


(3.46) 
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8. Solve the adjoint problem 

1 d{uX) divX) 


Pe dy^ 


+ 


dx 


+ 


dy 


+ 2{9-Y)6{y-yi) = 0 


with 


9. Compute the gradient 


a; = 1, A = 0 

y = 0, ^ = -2wPeie^-Y^) 
dy 

y = d, A = 0 
A(x, 0) 


J' 


Pe 


10. Compute the descent direction 

pk _ ^kpk—l jfk 

where 

for the first iteration 7 * = 0 
otherwise 7 = 


||J*=- 


11. Solve the sensitivity problem 


with 


d{uA9) d{vA9) _ 1 d^jAP) 

dx dy Pe dy"^ 


x = 0, A9 = 0 
y = l, A9 = 0 

y = 0, ^^^ = -ANu = P^ 

dy 


12. Compute the step size 

A9{9 - Y) ^Av - Vo) + ^iv) dx 




dx 


Jo - »;)+%) 

where mp indicates number of measurement planes and Wj is 
assigned to the data from the j-th plane of measurement. 


(3.47) 

(3.48) 

(3.49) 

(3.50) 

(3.51) 

(3.52) 

(3.53) 

(3.54) 

(3.55) 

(3.56) 

(3.57) 

(3.58) 

(3.59) 
the weight 


13. Compute the new local Nusselt number 

(_Nu)^+V= (-Nu)*^ - 


(3.60) 


14. Go to 5 and continue till convergence. 
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3.3 Steady Three-Dimensional Temperature Field 

The formulation of Section 3.2 is now presented for a steady three dimensional 
temperature over a flat plate. 

3.3.1 Assumptions 

1. Temperature field is three-dimensional, i.e. the plate temperature is varying 
with X and ^ coordinates 

2. The temperature field is steady. 

3. Boundary-layer approximations are applicable for the energy equation. 

4. The flow field continues to be two dimensional and the z component of 
velocity is zero. 

3.3.2 Direct Problem 

The direct problem solves the following energy equation in dimensionless form: 


d{ue) 

dx 

d{vd) _ 1 dH\ 

dy Pe \dy‘^ ) 

(3.61) 

The boundary conditions are: 

11 

O 

II 

o 

(3.62) 


y = 0, = -Nu 

dy 

(3.63) 


11 

11 

o 

(3.64) 



(3.65) 


. = r, ^=0 

(3.66) 


where, Nu is the local Nusselt number distribution. 
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3.3.3 Inverse Problem 

The inverse problem is posed in the form of minimization of the following func- 
tional: ^ 

J= f f [i^vj - Y^f + w{Bm - Ymf]dxdz (3.67) 

Jo Jo 

where Oiu and are the computed and measured temperatures at the wall, dm 
and Ym are the computed and measured temperatures in the interior and w is 
the weight assigned to the high quality data. 


3.3.4 Sensitivity Problem 


The sensitivity problem determines the the dependence of the temperature on 
the local Nusselt numbers. To formulate the sensitivity problem, a perturbation 
of is given to 9 and 9 is replaced by 0 + A^ in Equation 3.61 to yield 


d[u{9 + Ag)] g[u(g + Ag)] ^ 1 
dx dy Pe 


d'^{e + M) d^9 + A9) 


dy^ 




(3.68) 


Now, subtracting Equation 3.61 from Equation 3.68, we have the following sen- 
sitivity problem: 


d{uA9) d{vA9) _ J_ ( d‘^{A9) ^ g^(A^) >| 


dx 


dy 


Pe V 


+ 


dz‘^ ) 


(3.69) 


Similarly, in the boundary condition, Nu is replaced by Nu+ANu and 9 by 9-{- A9\ 
thus we have: 


a: = 0, 9 + A9^Q 

(3.70) 

0, = (Nu+ ANu) 

dy 

(3.71) 

y ~ d^ 9 A9 — 0 

(3.72) 

oz 

(3.73) 

a(_e + Ae) „ 

,=r,—^ = 0 

(3.74) 


Subtracting Equations 3.62, 3.63, 3.64, 3.65, 3.66 from Equations 3.70, 3.71, 
3.72, 3.73, 3.74, respectively, we get the following boundary conditions for the 
sensitivity problem: 

x = 0, A0 = O (3.75) 
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= 0, ®f)=-ANn 

dy 

(3.76) 

y = d, A9 = 0 

(3.77) 

z - n n 

(3.78) 

’ dz 

(3.79) 


3.3.5 Adjoint Problem 


Let, A be the Lagrange multiplier that is associated with the constraint. We can 
write the functional 

J = f r + w{dm - Ymf] dxdydz + 

Jo Jo Jo 

d{uB) d(v6) 


1 pd pT 



A 


0 Jo Jo 


Pe dz^ 


dx 


dy J 


dxdydz (3.80) 


In Equation 3.80, setting 6 + A0 in place of 6 and J + A J in place of J, we 


get 


J + AJ = f f f [{& + -yTS{y) + w{6 + A9 -Y)'^5(y - yi)]dxdydz + 
Jo Jo Jo 

dxdydz — 
dxdydz 


d pr 



0 JQ 
X pd pr 


A 

Pe 


d\9 + j\e) a^(g + Ag) 



d{u{9 + A0)} d{v{6 + A6)} 

_ +. _ 


(3.81) 


'0 JO JO 

where, 6{y — yo) is the Dirac delta function centered around yo 

Subtracting Equation 3.80 from Equation 3.81 and neglecting the terms in- 
volving the squares of A0, we get 

r»l pd pr 


AJ - 2 [ [ [ A9[w{9 -Y)6{y) + {9 - Y)5iy- yi)]dydxdz + 
Jo Jo Jo 

dxdydz — 


m *" A 

r*l pd pr 



Q Jo JQ 


d^A9) d^A9) 

dy^ dz^ 
d{u{A9)} d{v{A9)} 

dx dy 


dxdydz 


(3.82) 
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Using the boundary conditions from the sensitivity problem, we have 


) dv^ 


= A / -hr^dy 


d^jAd) 

9y2 


= A 


r<^ d{A6) dX 
Jy=a JO ^y 


r fdX fd^A9)^ \ 1 ''=^ 

diA») 8X 


y=d pd q2\ 


Therefore 


0 Jo Jo 


d{Ae) 


X —^ dxdydz 

dy2 


9y\„^ 


y=0 ^0 


5(A«) 


';(Ae)dy (3.83) 


dxdz + 


A9 dxdz + 


IQ JO i^y jy=o 

/•l pd pr a 1 \ 

L 1 1 


(3.84) 


Again 


• ^(A^) _ r 

X——dz - ^A 

= A^ 


d‘^{A9) 


d{A9) 


- A9- 


z=^r nr 


^ rax f d'^{A 9 ) 

\dz j dz^ ' 

d{A9) ax 


dz 1 dz 


2=0 Jo 9z az 
pr a2\ 

i„u/ 


(3.85) 


Therefore 




Jo Jo [dy 

r»l rd rT 


A9 dxdy + 


m r a 2 \ 

Ae^d.dyd. 


(3.86) 


/■‘A^h^drr 




9A r d{u{A9)} 
ax . ax 


[An(AC3- 

/•l Q\ 

[AuA^l^^i - J^ u-^{A9)dx 


(3.87) 
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Therefore 


I I I = j‘ l\xuAe]„,dydz 


-I rd />r Q\ 

u-^{A6)dxdydz (3.88) 



0 JQ JQ 


Further 




/O 


dy 


I 




dy J \dy J dy 

= - I" 


y=d 

2/=0 


(3.89) 


Therefore 


I! [ r =-/'/' [ v^mdxdydz (3.90) 


1 pd pr 


'0 JQ JQ dy 

Combining Equations 3.82, 3.84, 3.86, 3.88, 3.90, we have 


1 pd pr 


AJ = 



'0 Jo JO 
''I pd pr 


1 fd'^X d‘^X\ dX dX 
Pe \dy^ dz^ J'^^dx^^dyj 


Addxdydz + 

2 f [ [ [w{9 - Y)5{y) + {6- Y)S{y - yi)] ABdxdydz + 
Jo Jo Jo 




dX 


dxdz + 


i_ [' f 

Ps Jo Jo 


X 


d{Ae) 
dy J 


-) j/=d 


dxdz — 


y=0 


A9 


dz 


dxdy 


(3.91) 


z-0 


dyjy^o 

pd pr pi pd 1 

ll[Xudsei,,dydz-ll- 

For AJ ->■ 0, the following adjoint equation and boundary conditions can be 
deduced, 

(3.93) 

(3.94) 

(3.95) 

(3.96) 


1 

Pe 


a: = 1, A = 0 

y = 0, ~ = -2wPe{9^-Y^) 

dy 

y = d, X = 0 
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(3.97) 

Hence 



A J = ^ ^ [A(- AJVu)]^^o dxdz 

(3.98) 

From the definition of a gradient, we may write. 

AJ = r [ J'{-Am)dxdz 

Jo Jo 

(3.99) 


From Equations 3.98 and 3.99, we have the following expression for the gra- 
dient of the functional; 


J’ = lim 


AJ A(a;, 0,2) 


aNu-^-o ANu 


Pe 


(3.100) 


3.3.6 Determination of Search Step Size. 

If the value of the functional is J* after the fc-th iteration and it becomes 
after temperatures are updated through a step size of ,5*, the condition for the 
optimum step size will be given by 




0 


(3.101) 


Recalling the expression of functional 

pd pr 


Jk ^ f" f’' r _ Y)\y)] dxdydz (3.102) 

Jo Jo Jo 


we get 


mm„j‘(-AiV«)‘ = mm„ f" f {wle{[- AN uf - )3‘P‘} - YfKv - Vi) 

JQ Jo Jo 

+[e{{-ANu)’‘-/3'^P'^}-Yf6{y)}dxdydz (3.103) 


Linearizing Equation 3.103 by Taylor’s series expansion 

pd pT 


== w ^ ^ {e[{-ANuf - ^*A0(P'')] - YY6{y - yi)dxdydi 

Y t t - ^’^A9{P^)]-Y}H[y)dxdydz 

Jo Jo Jo 

m {(^ - Y) - ^’'AdY{w5{y - yi) + 5{y)}dxdydz 

■ ' ' ' ■ 


(3.104' 



3.3 Steady Three-Dimensional Temperature Field 


42 


Therefore, 

a/3* 

Using Equations 3.101 and 3.105, we have the following equation for the op- 
timum step size, 

Qk.^ fo fo^ fo - yi) + ^o^(p)}dxdydz ^ 

fo lo - yi) + 

3.3.7 Complete Inverse Algorithm. 

The complete inverse algorithm for calculation the Nusselt number distribution 
over flat surface is summarized below: 

1. Solve the continuity and momentum equations to get the two dimensional 
velocity field. 

2. Filter the imperfect plate temperature data using Gram orthogonal poly- 
nomial filter. 


m r 

A0(e - y - j3'^A&){S(y - yi) -+■ wd(y)}dxdydz (3.105) 

. 


3. Solve the three dimensional energy equation using the following boundary 
conditions: 


o 

II 

o' 

11 

(3.107) 

y = 0, 6 = 6^ 

(3.108) 

o 

II 

II 

(3.109) 


(3.110) 

n 

(3.111) 


4. Get the values of local Nusselt numbers. Take an initial guess of Nusselt 
number based on the distribution obtained. Initialize J = 0 

5. Solve the three dimensional energy equation using the following boundary 

conditions (with Nusselt number from step 4): 


a: = 0, 0 = 0 


(3.112) 
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de 

y = 0, — = -Nu 
dy 

(3.113) 

o 

II 

11 

(3.114) 

o 

11 

o 

11 

(3.115) 

de 


(3.116) 


Find the temperatures at the plate as well as the plane of high quality data. 


6. Compute the functional 

rd /*r 


m r 

{{0^ - + w{6m - Ym)^}dxdydz (3.117) 

7. Check the stopping criterion 


where, e is a predetermined real number, very close to zero. 
If it is satisfied, stop iterations. Otherwise 


8. Solve the adjoint problem 


1 d(uX) d(vX) , , . 

Sx + Sy 

(3.119) 

with 


X = 1, A = 0 

(3.120) 

o \ 

y = 0, — = -2wP&{dy, - Yy,) 

(3.121) 

y = d, A = 0 

(3.122) 

^ dX 

(3.123) 


(3.124) 

P. Compute the gradient 

rf A(x,0,2:) 

Pe 

(3.125) 
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10. Compute the descent direction 


pk ^ ^kpk-l ^ J,k 

(3.126) 

where 


for the first iteration 7 * = 0 

(3.127) 

otherwise 7 = — - 

(3.128) 

11 . Solve the sensitivity problem 


d{u/\9) d{vA6) _ 1 ra2(A0) d\Ady 
dx ' dy Pe [ dy'^ ' 

(3.129) 

with 


x = 0 , = 0 

(3.130) 

^ = 0, = ANu = 

oy 

(3.131) 

y = d, A9 = 0 

(3.132) 


(3.133) 

Sm n 
" = az =“ 

(3.134) 

12 . Compute the step size 


^ fo So [^T=i ^3^ (y ~ 

^ - /; /;(Ae)2 [Er^i - ») + «(!/)] 

(3.135) 


where mp indicates number of measurement planes and Wj is the weight 
assigned to the data from the j-th plane of measurement. 


13. Compute the new local Nusselt number 

(_Nu)*+^ = (-Nu)^ - (3.136) 

14. Go to 5 and continue till convergence. 
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3.4 Unsteady Three-Dimensional Temperature 
Field 

The inverse estimation of the instantaneous Nusselt number distribution over a 
flat surface when the thermal field is three dimensional is presented in this section. 
The flow field is taken to be steady and two dimensional. 


3.4.1 Direct Problem 

The direct problem solves the following energy equation: 

W d{ue) djve) _ J_ ^ 

dt'^ dx dy Pe **" dy'^ dz^j 
The initial and boundary conditions are: 



(3.137) 

(3.138) 

(3.139) 

(3.140) 

(3.141) 

( 3 . 142 ) ; 

(3.143) 

(3.144) 


3.4.2 Inverse Problem 

The inverse problem is posed in the form of minimization of the following func- 
tional: t 1 /T 

J.^{\ \\{Q^-Y^f^w{Qyj^~Yr^^\dxdzdt (3.145) 


/o ao ao 


[(^,0 - -~ Ymf]dxdzdt 


where and are the computed and measured temperatures at the wall, 0^ 
and Ym ^re the computed and measured temperatures in the interior and w is 

the weight assigned to the high quality data. 
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3.4.3 Sensitivity Problem 


The sensitivity problem determines the the dependence of the temperature on 
the local Nusselt numbers. To formulate the sensitivity problem, a perturbation 
of Ad is given to 9 and 9 is replaced hy 9 + A9 in Equation 3.137 to yield 


d{9 + A9) d[ui9 + A0)] d[v{9 + A0)] _ 

dt dx dy 

i_ + + Ag) d'^{9 + A9y 

Pe _ dx"^ dy"^ dz^ 


(3.146) 


Now, subtracting equation(3.137) from equation(3.146), we have the following 
sensitivity problem 


djAd) d{uA9) djvAe) _ J_ r a^(Ag) d^A9) d'^jAd) ' 

dt dx dy Pe _ dx'^ dz^ 


(3.147) 


Similarly, in the initial and boundary condition, Nu is replaced by Nu+ANu and 
9 hy 9 A9] thus we have: 


f = 0, 0 + A0 = O 

(3.148) 

X = 0, ^ -p A0 = 0 

(3.149) 

dy 

(3.150) 

= -(Nu + ANu) 
ay 

(3.151) 

y = d, 9 + A9 = 0 

(3.152) 

. = 0. + 
az 

(3.153) 

az 

(3.154) 


Subtracting Equations (3.148-3.154) from Equations (3.155-3.160) respectively, 
we get the following boundary conditions for the sensitivity problem; 
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, = 0,^ = -ANu 

oy 

y = d, A0 = 0 


0 




dz 


(3.158) 

(3.159) 

(3.160) 

(3.161) 


3.4.4 Adjoint Problem 


Let, A be the Lagrange multiplier that is associated with the constraint. We can 
write the functional 


J 


r [ P r X 
to Jo Jo Jo 

ntf pi pd pT 


1 (dH d'^d d^9\ d{ue) d{v9) 89 

+ + 


[Pe \dx^ dy^ 9^2 


dx 


dy dtj 


/ f f f [i^w -y-w)'^ + w {9m- ymY]dxdydzdt 
Jo Jo Jq Jo 


dxdydzdt 

(3.162)' 


In Equation 3.162, setting 9 + in place of 9 and J + A J in place of J, we 
have 


J+AJ = 

rtf pi pd pr 


f f f f [{9 + A9 — y)^5{y) + w{9 + A9 — yy6{y — yi)]dxdydzdt 
Jo Jo Jo Jo 

rtf rl rd rr \ i A ^ A Q\ 5;i2/^/D i A /3M 

JO JO Jo t/o 


+ 


1 pd pr ^ 

Pe 

4f pi pd pr 

A 


’ d^{9 + A9) d'^{9 + A9) d^{9 + A9) 

dx^ dy^ 8z^ 


dxdydzdt 



to Jo Jo Jo 


d{u{9-^Ad)) d{v{9 + A9)} d{9 + A9) 

dx dy dt 


dxdydzdt 

(3.163) 


where, 5{y — yo) is the Dirac delta function centered around y^ 

Subtracting Equation 3.162 from Equation 3.163 and neglecting the terms 
involving the squares of A^, we get 

AJ = 2 f ^ [ [ f A9 [w{6 - y)6{y) + {9 - y)5{y - yi)] dydxdzdt 
JO JO Jq Jo 

pf p r A 
Jo Jo 


: + 



to Jo Jo Jo 1^® 

Jtf nl rd rr ^ 

Jo Jo Jo Jo 


' 8^A9) d^A9) d^{A9) 

dx^ 9y2 9^2 


dxdydzdt — 


9{«(A^)} , 9{u(A0)} , d{A&) 

— ^ 1 . 1 


dx 


dy 


+ 


dt 


dxdydzdt 


(3.164) 
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Using the boundary conditions from the sensitivity problem, we have 

/ 

d{Ae) dx 

Jo 


,9"(A9) 

. = 


d{A9) 


x=0 


dx 

d{Ae) 

dx 


J a:=0 
x=l 


x=0 


dx dx 


Ae 


dx 


1 X=1 


+ 


d'^X 


J a;=0 


r a(A9)i 



dx 

x=0 

dx_ 


+ / 

J x=l 0 


0 dx"^ 

^ d‘^X 
dx^ 


{Ae)dx 


{A6)dx 


(3.165) 


Therefore 

rtf rl nd f*r 



'0 Jo Jo Jo 


X -^P dxdydzdt 
dx^ 


if fd 



X 


'0 JO JO L 
rtf pd pr 

dx 


d{Ae) 

dx 


-L 

m r - 

Ja;=l 

ptf pi pd pr op,\ 

I ill 


dydzdt — 


Ae 


x=0 

dydzdt + 


(3.166) 


Again 




Xl^iy- 


dy^ 

d{Aey ^=^ 

dy 


dx fd^{Ae) ^ 


JW*)* 


A 

-=o Jo 


dy J dy^ 

^ d{Ae) dX 


Jy=0 


2/=0 

y:zzd 


d{Ae) 

9y Jj;=o 
djAP) 
dy 


y- 

y-d 


Ae 


+ 


yz=Q 


dy dy 
dy 


dy 


y=d nd fs2\ 
pd «2 \ 


(3.167) 
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Therefore 




m r r 


A -— — - dxdzdt + 


dy J^=o 


tf rl rr 



0 Vo JO L 


dy 


Ae 


dxdzdt + 


Jy=o 


ptj pi pd pr o2\ 

L III 


(3.168) 


Again 


= 

az^ 




J z=0 


A 


d{A6) 


dz 


d{Ae) dx 


Ae 


dX 

dz 


Jz=0 Jo 


I 


z=o JO dz dz 

z=r rr q2^ 

+ ' 


(3.169) 


Therefore 

pf p p p d^{Ae) 

A" 



'0 Vo Vo Vo 


dz^ 


dxdydzdt 


m d 


dy 


Ae 


dxdydt + 


J z=Q 


.1 M rr 


I'lhL 


(3.17D) 


Also 


I 


\aHA9)},, _ 


ax 


i 


d{u{Ae)} 


X I ■'^-p-^J- dx 
dx 


/{£/ 

/*! f)\ 

= iXu{Ae)rJ,-l 2iK(AS)& 

dX 

= [XuAei^^- u—iAe)dx 


dX f d{ u{Ae)} \ ' 
dx J 


dx 


X=1 


x=0 


(3.171) 


Therefore 


rtf ri I'd rr ^ d{u{Ae) }^^^^^^^^ ^ f f [AwA0]^^i dydxdt 
/o Jo Jo Jo >^0 


df fd r dX....,,,,. 

u^{A9)dxdydzdt 

0 Vo Vo VO 

(3.172) 
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And 


h dy 


X f 


dy 




= [Xv{Ae)]lJ - ~v{Ae)dy 


y=0 




(3.173) 


Therefore 


r‘‘ p !■“ r\ d{v(A})} 

'0 Jo Jo Jo 


rtf rl na nr 

dxdydzdt = — / / / v~{A6)dxdydzdt 

Jo Jo Jo Jo 

(3.174) 


Further 


I 


0 dt 




t—tf 

t=0 


ptj pj\ 

[AA^lSS^ - Ae^dt 


rtf 


(3.175) 


Therefore 


^tf ni nd r ^ rf r [xAe\^. dxdydzdt 

Jo Jo Jo ^ Jo Jo Jo Jo ’ 


tf nl nd 


'0 v^O -'0 ^0 


[/ [ [ r Ad^dxdydzdt (3.176) 

Jo Jo Jo Jo 

Combining Equations 3.164 to 3.176, we have . 

AJ= ■ 

ntf nl nd nr f / g2 ^ 02 X d^X\ , dA , SX , dAl ... , , ' ^ ^ 

2 [ f [ [w{e-Y)5{y) + {9-Y)5{y~yi)]A9dxdydzdt- 

Jo Jo Jq Jo 


Hf rd 


If 


Pe 


m d ^ 

Pe [dy 

m r , 

[AAi9]^^t dxdydz 


A 


d{A9) 


dx 


dydzdt 


J a;=0 
1 z=tr 


ntf fd fr 

Pe 

rtf fd rr 


m r 


A9 


dx 


dydzdt 


X=1 


- z=tr f^tf pd pr 

A9 dxdydt- / / [Xu A9]^^^ dydzdt 

. 2=0 Jo Jo Jo 


(3.177) 
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For AJ 0, the following adjoint equation and boundary conditions can be 

deduced: 



-yi) = 0 (3.178) 

t = tf, A = 0 

(3.179) 

II 

O 

11 

o 

(3.180) 

X — 1, TT- + uAPe = 0 
ox 

(3.181) 

Q \ 

y = 0, — = -2u;Pe(^^ - Y^) 
ay 

(3.182) 

O 

11 

11 

(3.183) 


(3.184) 


(3.185) 

Hence ^ 

= [Xi^ANu)]^^,dxdzdt 

(3.186) 

From the definition of gradient, we may write, 


AJ = r / ' r J'(-ANn)dxdzdt 

Jo Jo Jo 

(3.187) 

From Equations 3.186 and 3.187, we have the following expression for the 

gradient of the functional: 


, AJ X{t,x,0,z) 

^ ' aNu-o A(Nu) " Pe 

(3.188) 


3.4.5 Determination of Search Step Size. 

If the value of the functional is J'' after /c-th iteration and it becomes after 
temperatures are updated through a step size of the condition for the optimum 
step size will be given by 




Recalling the expression of functional 

ptf rl nd rr 

Iq Jo Jo Jo 




[{6 - Yf5{y - yi) + w{6 ~ T)<5(t/)] dxdydzdt (3.190) 
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Therefore, 


min^ r f' f" - YfSiy - yi) 

Jo Jo Jo Jo 

+w[d{{~ANu)’^ - - Yf5{y)}dxdydzdt 

(3.191) 


linearizing by Taylor series expansion 


jfc+i[(_AiVw)^+^] = 

/ V / F [^[(-AiVu)* - P'^AeiP^)] - Yf5{y - y{)dxdydzdt 
Jo Jo Jo Jo 


+W [^ [ f [ [9[{-ANu)^ - /3'^A9{P’‘)] - Y? 6{y)dxdydzdt 
Jo Jo Jo Jo 

f f f f [{^ -y) - ^'‘^0V{^{y-yi) + w5{y)]dxdydzdt 

Jo Jo Jo Jo 


(3.192) 


Therefore, 




otf p\ pd pT 

= 2 / / / / A9{9—Y A9){5{y—yi)-\-w5{y)}dxdydzdt (3.193) 

Jq Jo Jo Jo 


Using Equations 3.189 and 3.193, we have the following equation for the op- 
timum step size, 

Rk ^ fo^ Jo fo So ~ dxdydzdt 

/o^ /oVoVo’'(A^)M<5(y - yi) + 


3.4.6 Complete Inverse Algorithm. 

The complete inverse algorithm for calculating the Nusselt number distribution 
over a ribbed surface is summarized below: 

1. Solve the continuity and momentum equations to get the velocity field. . 

2. Smooth imperfect plate temperature data using Gram orthogonal polyno- 
mial filter. 

3. Solve energy equation using the following initial and boundary conditions: 

i = 0,0 = 0 (3.195) 
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a; = 0, 9 = 0 

(3.196) 


(3.197) 

II 

o 

II 

(3.198) 

y = d, 9 = 0 

(3.199) 

^89 

(3.200) 


(3.201) 


4. Get the values of local Nusselt numbers. Take an initial guess of Nusselt 
number to be same as that found at the end point of the domain. Initialize 
J = 0 

5. Solve energy equation using the following initial and boundary conditions: 


t = 0, 9 = 0 

(3.202) 

a: = 0, 9 = 0 

(3.203) 


(3.204) 


(3.205) 

tl 

"" < 

II 

O 

(3.206) 

. 89 ^ 

(3.207) 

■ 99 ^ 

(3.208) 


Find out the temperatures at the plate and at the plane of measurement of 
redundant data. 

6. Compute functional, 

J = f ^ f f f {{e^ - + wiSm - ymf jdxdydzdt (3.209) 

Jo Jq Jo Jo 

7. Check the stopping criterion, 

(3.210)- 

where, e is a predetermined real number, very close to zero. 

If it is satisfied, stop iterations. Otherwise: 
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8. Solve adjoint problem 

1 . 4. ^ ^ 


Pe dz^ J dx dy dt 

with 


+ — + 2(0 - F)(5(j/ - yi) = 0 



(3.211) 

t = tf, A = 0 

(3.212) 

rr = 0, A = 0 

(3.213) 

. dA 

X = 1, 7: — h uAPe = 0 
ox 

(3.214) 

0, ^ = -2tePe(0^u - y^) 

(3.215) 

0 

II 

^ ; 

II 

SJs 

(3.216) 

0 

II 

:C 5 |cc> 

o' 

11 

(3.217) 

„ 

. = r. jj = 0 

(3.218) 

7/ _ A(i, 5 ;j 0 , 2?) 

Pe 

(3.219) 


9. Compute gradient 


10. Compute the descent direction 

pk^^kpk~i_^j,k 

where 

for the first iteration 7*^ — 0 

r\T h Ci-rxxn o/i — -J ! L 


otherwise 7 

11. Solve the sensitivity problem 


i4‘-'ir 


(3.220) 

(3.221) 

(3.222) 


a(A«) , a(«A9) , a(«Ae) _ 1 r8=(Ae) s=(Ae) s^(a9)1 


t = 0, A0 = 0 
X = 0, A0 = 0 
a;-l 

Jj — ij Q — 


(3.224) 

(3.225) 

(3.226) 
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, = ANu = P* 

dy 

(3.227) 

y = d, A9 = 0 

(3.228) 

. = 0, 

oz 

(3.229) 

z = r, ^ = 0 

OZ 

(3.230) 


12. Compute the step size 


/o' /o s; - y) [E”Si - yj) + Hv)[ 

dxdzdt 

/„" /o [Eri - Vi) + Hv) 

dxdzdt 


(3.231) 


where rap indicates number of measurement planes and Wj is the weight 
assigned to the data from the j-th plane of measurement. 


13. Compute the new local Nusselt number 

(_Nu)''+^ = (-Nu)* - (3.232) 


14. Go to 5 and continue till convergence. 



Chapter 4 

Flow and Heat Transfer over a 
Surface with a Rib 


The mathematical details of the inverse technique based on iterative regulariza- 
tion is explained in this chapter for flow and heat transfer over a ribbed plate. 
The basic algorithm is similar to Chapter 3. However, due to presence of rib, the 
boundary layer approximations are no longer valid here. The direct, adjoint and 
sensitivity probleips need special treatment at the rib location. 

The Figure 4.1 shows the single rib configuration that has been analyzed. 
The rib cross-section is rectangular and the cross-section does not change in 
the z direction. The fiow field is taken to be steady two-dimensional, while 
unsteady two and three-dimensional temperature fields have been considered. 
The following assumptions have been employed throughout: 

1. The fluid is Newtonian and incompressible. 

2. The fluid properties are constant with respect to temperature. 

3. The flow is laminar. 

4. The velocity field is steady and two-dimensional. 

5. The fluid inlet temperature is spatially uniform and temporally constant. 

6. Viscous dissipation and mechanical work terms are negligibly small. 

7. Forced convection is predominant; hence the momentum and energy equa- 
tions are not coupled. 
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4 K 



Figure 4.1: Plate with Single Rib 

8. The flat plate temperature is temporally a constant. 

9. Plate temperature data is available with some scatter whose characteris- 
tics are known; in addition, limited but high quality temperature data is , 
available within the flow field. 

For steady surface temperature case the algorithm consists of the following 
basic steps: 

1. Smooth the given temperature data using Gram orthogonal polynomial 
filter through the seven point moving window averaging scheme. 

2. Solve the continuity and momentum equations to find the velocity field. • 

3. Take an initial guess of local Nusselt number distribution over the plate. 

4. Compute the temperature field using a guessed local Nusselt number as the 
boundary condition at the plate. 


5. Compute the optimization functional. 
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6. If the value of the functional conforms to the stopping criteria, stop itera- . 
tions and present the results. 

7. Otherwise, solve the adjoint problem and find the gradient of the functional. 

8. Find the search direction. 

9. Solve the sensitivity problem and find the search step size. 

10. Update the local Nusselt number values over the plate surface. 

11. Return to 4 and continue. 


4.1 Velocity Computation 


The two-dimensional velocity field is computed by solving continuity and Navier- 
Stokes equations through the vorticity-stream function approach. The governing 
equations are: 

( 4 . 1 ) 

( 4 . 2 ) 

( 4 . 3 ) 

( 4 . 4 ) 

( 4 . 5 ) 

( 4 . 6 ) 

( 4 . 7 ) 

( 4 . 8 ) 

( 4 . 9 ) 
( 4 , 10 ) 


doj d(uuj) d(vaj) _ 1 / d'^cj 

dt dx dy Re \ dx"^ dy'^ 

d'lj) 

d'ljj 

The boundary conditions are: 

0, w = 0, tp = y, u — l, V = 0 


X 


_ 7 _n ^-n — = (] 
^ dx ~~ aa:2 '' ^ 


dx 


dx 


0 


2/,= 0, w = -VV) = n = 0, u = 0 
y = h, 00 = ~V'^tp, ‘tp ~ y, u = l 
X = li, and y < drii w ip = 0, u = 0, u = 0 

X = -b ?ri, and y < dri. w = V* = 0> “ = = 0 

x<li + lri, Siud X > h, andy = drii w = ^ = 0, u = 0, v = 0 (4.11) 

The rib is treated by the following condition: 

x <li + lrii and r > ^ 1 , and y < dri) ^ — 0) V' = 0) oj = 0 (4.12) 
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4.2 Unsteady Two Dimensional Temperature Field 

The inverse formulation for a two dimensional unsteady temperature field is de- 
scribed here. The flow field is considered steady. 

4.2.1 Assumptions 

1. Plate temperature is steady. 

2. Plate temperature is varying with length only. 

3. Temperature field is developing. 


4.2.2 Direct Problem 


The direct problem solves the following energy equation: 

dt dx "** dy Pe J 

The initial and boundary conditions are: 


(4.13) 


t = 0, 9 = 0 

x = 0, ^ = 0 
X = 1, 9 = 0 

m 

dy 

y = d, 9 = 0 


y = o, 


-Nu 


(4.14) 

(4.15) 

(4.16) 

(4.17) 
{4.18) 


where, Nu represents local Nusselt number values (initial guess or updated values 
after any iteration). The rib is treated by the following condition. 


For X <li + Iru ^ ^ ^i> y ^ 
ag 1 g, fd'^9 cf^9 

89 


For X = li or x = li-\- hu ^.nd y < dri, 

89 

For X > h or x < h + Iri) ^.nd y = dru ^ 


solid 


= kf 


8x 
= kf^ 

solid 


fluid 

fluid 


(4.19) 

(4.20) 

(4.21) 




4.2 Unsteady Two Dimensional Temperature Field 


60 


4.2.3 Inverse Problem 


The inverse problem is posed in the form of minimization of the following func- 
tional; t i 

~ J J - yw)"^ + w{6rn - Yrn)'^]dxdt (4.22) 

where and Yyj are the computed and measured temperatures at the wall, 9m 
and Ym are the computed and measured temperatures in the interior and w is 
the weight assigned to the high quality data. 

4.2.4 Sensitivity Problem 

The sensitivity problem determines the the dependence of the temperature on 

the local Nusselt numbers. To formulate the sensitivity problem, a perturbation 

of A0 is given to 9 and 9 is replaced by 0 -f- A0 in Equation 4.13 to yield 

d{9 + A9) d[u{9 + A9)] d[v{9 + A0)] 

H b ' _ 


dt 

1 

Pe 


dx dy 

d^{9 + A9)^d‘^{9 + A9) 


dx^ ' J 

Now, subtracting Equation 4.13 from Equation 4.23, we have the following sen- 
sitivity problem, 


d{A9) d{uA9) d{vA9) ^ 1 

dt dx dy Pe 


d\A9) d^iA9) 


(4.24) 


dx'^ dy"^ J 

Similarly, in the initial and boundary condition, Nu is replaced by Nu+ANu and 
6* by ^ + A9\ thus we have: 

t = 0, ^ + A^ = 0 (4.25) 

X — 0, 9 -j- A9 = 0 (4.26) 

(4.27) 


dy 


y 


0 = _(Nu + ANu) 

dy 


■ (4.28) 

y = d, 9 + A9 = 0 (4.29) 

Now, subtracting Equations ( 4 . 14 - 4 . 18 ) from Equations (4.25-4.29) respectively, 

we have the following boundary conditions for the sensitivity problem 

t = 0, A9 = 0 (4.30) 
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a: = 0, = 0 


a; = 1, — = 0 

X 


y = o, 


d(A0) 


y = d, A(9 = 0 

The rib is treated by the following condition. 


(4.31) 

(4.32) 

(4.33) 

(4.34) 


For X < li + Iri, and x > li, and y < dri 
d(Ae) ^ f d^jAe) d^Ae) \ 
dt Pea/ \ dx^ dy^ J 


(4.35) 


For x = li ov x = h + lri, and y < d^i, ks d(A^ (4.36) 

dx solid fluid 


For X >li or X < h + lri, and y = dri, ks = kf I (4.37) 

9y solid I- - 


d{A9) 

fluid 


4.2.5 Adjoint Problem 


Let, A be the Lagrange multiplier that is associated with the constraint. We can 
write the functional 

J = f' J f\(e„-Y„f + w(e„-Y„)^]dxdydt + 

r' /■« ri fdH , s^e , dH\ a{«e) d{ve) 


If 


^ [Pe ~di ~ ^ 

(4.38) 


In Equation 4.38, setting 6> + A0 in place of 9 and J + A J in place of J, we 


J+AJ = 

ptf pi pd 




'0 v/O ^0 
ptf pi p 


Pe dx^ 


f f ^ ^ ^ ^ ^ H ^ ^ ^ I dxdydt 

Jo Jo [ Sx dy at J 


C4,39) 
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Subtracting Equation 4.38 from equation 4.39 and neglecting the terms in- 
volving squares of A6, we have 

ptf pi pd pr 


AJ 


2 1 III Ad[w{e~Y)5{y) + {9-Y)S{y-yi)]dydxdzdt + 
Jo Jo Jo Jo 

rtf pi pd pr ^ 


pif pL pa pr 

Jo Jo Jo Jo 


Pe 


d^{Ae) (AO) 


dx"^ 




dxdydt 


m d 

_ 


A 


5{u(A^)} d{v{Ae)} d{A&) 


dx 


- 1 - 


dy 


+ 


dt 


dxdydt 


(4.40) 


Using the boundary conditions from the sensitivity problem, we get 

d{Ae)Y=^ f^d{Ae)dx 


X 


am ,, = 


dx"^ 


dx 

d{Ae) 
dx 

d{Ae) 


x—0 

x^l 


2=0 


Jo 


A9 


dx dx^"" 

dAr=' . f^d'^X 


X 


dx 


J ic=0 


dx 

A9 


X’ 

dx 




(4.41) 


Therefore 


tf p Y d^{A9) 


fl 


X- 


dx"^ 


-dxdydt 


= _ r [a 

Jo Jo . 


d{A9) 

dx 


dydt ■ 


x^O 


n d - 

_ 


dx 


A9 


dydt -I- 


X=l 


j A4. 

A9 -^dxdydt 

dx 


(4.42) 
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Again 


, d^{Ae) , 

A--— A - dy 


dy^ 


Therefore 


5(^0) -jy=d 





d{Ae) dx 


Jy=Q 


dy dy 

-iV=d i.dQ2)^ 


‘^ d^x 

dy^ 


{Ae)dy 


(4.43) 




m r - 

- 


A 


y=d 


d{Ae) 

% J3;=0 


dxdt + 


/■'/[ 


dy 


Ae 


dxdt + 


2/=0 


m d q 2 \ 


(4.44) 


Also 


/' 


X 


d{u{A9)} 

dx 


dx = 


/'I 

= [At*(A«)]'i - ^n(Ae)dx 

/I 

= [AuA9]„, - j( «^(Ae)<fa 


1 x=i 
x=0 


(4.45) 


Therefore 


Jo JO JO '^o 


r I' j\^{Ae)dxdydt 


d\. 


(4.46) 
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Again 


= 


/o 


dy 


A 


d{n(Ag)} 

dy 


dy - 


dX fd{v{A9)} 

dyj 


dy 


dy 


^y=zd 

y=0 


fd 

= [At^(A9)lg - 


Therefore 




Further 


I 


dy 

d(A0) 

0 


(4.47) 


(4.48) 


dt 


A j ^^^dt 


J 

= [AA9]£‘' - r AS 

/*t/ 


a(A0) 
dt J dt 
*/ . 

dt 


A6dt 
dt 


t^tf 

J t=o 


(4.49) 


Therefore 


l‘ j\?^dxdydt = j'' l‘ j\xAe],,,^ dxdydt - 


'tf M nd 


m d 

AO-^dxdydt 
dt 


(4.50) 


Now combining Equations 4.40 to 4.50, we have 
AJ = 


■*/ /-n 1 fd^X d^X 




Pe Vdx^ dy^ 


dX dX dX 


AOdxdydt + 


dydt — 


2 t' t f [■m{e-Y)S{y) + {9- Y)5{y-yi)]Aedxdydt 

dydt- f fA\Aef 

h Jo Pe L 3x Jo Jo Pe L 9x,^, 

f‘' j\xuAe]^^, dydt- [XAe]^^,/xdy (4.51) 

For A J -4 0, the following adjoint equation and final boundary conditions can 
be deduced 

J_ f A + A A ?JAl + + ^ + 2(« - y)i(!/ - »i) = 0 (4.52) 

Pe Vax2 ^ dyV dx dy dt 
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t = tf, A = 0 

x = 0, A = 0 

X = 1, -^ + tiAPe = 0 
ox 

y = 0, — = -2toPe(^^ - Y^) 
y = d, A = 0 

The rib is treated by the following condition. 

For X < li + Iri, and x > li, and y < dri 
dt Pea/ \dx‘^ dy"^ 


(4.53) 

(4.54) 

(4.55) 

(4.56) 

(4.57) 


Since 


For x — h or a: = + Iri, and y < dri, K 
For x>li ox X <li+ Irii and y = dr\, k 

1 r’’ 


dx 

dy 


solid 

solid 


^kf 


= kf 


dx 

dy 


fluid 

fluid 


AJ = f' f (a(-a;v«)1„„, dxdt 

and the definition of gradient 

ftj pi 


A J — f f J'{—ANn)dxdt 
Jo Jo 


(4.58) 

(4.59) 

(4.60) 

(4.61) 

(4.62) 


From Equations 4.61 and 4.62, we have the following expression for the gra- 
dient of the functional: 


P- ■ AJ A(t,a:,0) 
"aNu->oANu Pe 


(4.63) 


4.2.6 Determination of Search Step Size. 

If the value of the functional is J*' after fc-th iteration and it becomes after 
temperatures are updated through a step size of 0^, the condition for the optimum 
step size will be given by 


d0^ 


= 0 


(4.64) 


Recalling the expression of functional 

rtf pi pd 


[{e ~ Yf5{y - yi) + w{e - y)5(y)] dxdydt (4.65) 

Jo Jo Jo 
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Therefore 


mine r ( [ {[^{(-AiVu)"' - - Y^y - y{) 

Jo Jo Jo 

+w[e{{-ANu)’^ - - Yf6iy)}dxdydt 

(4.66) 


linearizing by Taylor series expansion 


jfc+i[(_AiVu)*+^] = 

I' fo Jo h^[(-AiVu)'= - /3'‘A6{P>^)] - Yf5{y ~ y{)dxdydt 
~ lo lo 'Io + 


(4.67) 


Therefore 

dJ’^+^ 

dp’^ 


= 2 [' f [‘^ A6{e -Y- l3>=A9){w.5iy - Vi) + 5{y)}dxdydt (4.68) 
Jo Jo Jo 

Using Equations 4.68 and 4.69, we have the following equation for the opti- 
mum step size 


13^ = 


lo^ lo fo M(^/ - ^i) + ^(^)] dxdydt 

fo^ fo foi^^)^ H(y - yi) + <^(v)] dxdydt 


(4.69) 


4.2.7 Complete Inverse Algorithm. 

The complete inverse algorithm for calculating the Nusselt number distribution 
over a ribbed surface is summarized below: 

1. Solve the continuity and momentum equations to get the two dimensional 
velocity field. 

2. Filter the imperfect temperature data using Gram orthogonal polynomial 

filter, . 

3. Solve the three dimensional energy equation using the following initial and 
boundary conditions: 

t = 0, 0 = 0 


(4.70)- 
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X = 0, 9 = 0 


y = 0, 9 = 9 ^ 

y = d, 9 = 0 


The rib is treated by the following condition. 


X < li + Iri, and X > li, and y < dri, u = 0, -y = 0 


(4.71) 

(4.72) 

(4.73) 

(4.74) 

(4.75) 


4. Get the values of local Nusselt numbers. Take an initial guess local Nusselt 
number. Initialize J = 0 

5. Solve the two dimensional energy equation using the following initial and 
boundary conditions (with Nusselt number from step 4): 


t = 0, 6 = 0 
x = 0, 9 = 0 


89 

, = 0, jj = -Nu 

y = d, 9 = 0 


The rib is treated by the following condition. 


(4.76) 

(4.77) 

(4.78) 

(4.79) 

(4.80) 


x<li+lri, Si.ndx>h, andy<dri, U = 0, v = 0 (4.81) 

Find the temperatures at the plate as well as the plane of high quality data. 


6. Compute the functional 

r [‘ fU (K - y.)" + - Y„f}dxdvdt (4.82) 

Jo Jo Jo 

7. Check the stopping criterion, 

|j*_ j*-i| <e (4.83) 

where, e is a predetermined real number, very close to zero. 

If it is satisfied, stop iterations. Otherwise 
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8 . Solve the adjoint problem 

dx 




Pe \ dx"^ dy^ 
with 


dy dt 


t = tf, A = 0 

X = 0, A = 0 

. 9X „ 

X = 1, — + wAPe = 0 
ox 

y = 0, ~ = -2wPe{dio - 
y = d, A = 0 

The rib is treated by the following condition. 

X < li + Iri, and X > li, and y < dri, u = 0, v = 0 


9. Compute the gradient 


J' = 


X{t, X, 0 ) 
Pe 


10. Compute the descent direction 

pk _ ^kpk-l _|_ j/k 

where 

for the first iteration 7 ^ = 0 
. otherwise 7 *^ = 


_ ^jk _ jk~i^jk^ 




- 1||2 


11 . Solve the sensitivity problem 

a(Ag) djuAe) 9(uAg) ^ _ 1 _ 

dy Pe 


d\Ae) d^AO) 


dx^ 


dy^ j 


with 


X 


t = Q, A0 = O 
X = 0, A9 = 0 


y = 0 , 


d{A9) 


dy 

y — d} A6 = 0 


ANu = 


= 0 (4.84) 

(4.85) 

(4.86) 

(4.87) 

(4.88) 

(4.89) 

(4.90) 

(4.91) 

(4.92) 

(4.93) 

(4.94) 

(4.95) 

(4.96) 

(4.97) 

(4.98) 

(4.99) 
(4.100) 
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12. Compute the step size 


Ae(e - 

5 ^) 

Y^Tii^jHy - yj) + Hv) 

dxdt 

[> 

to 

- Vj) + ^{y) dxdzdt 


(4.101) 


where mp indicates number of measurement planes and Wj is the weight 
assigned to the data from the j-th plane of measurement. 


13. Compute the new local Nusselt number 

(— = (-Nu)*’ — 

14. Go to 5 and continue till convergence. 


(4.102) 


4.3 Unsteady Three Dimensional Temperature 
Field 

The inverse estimation of the instantaneous Nusselt number distribution over a 
flat surface when the thermal field is three dimensional is presented in this section. 
The flow field is taken to be steady and two dimensional. 


4.3.1 Direct Problem 

The direct problem solves the following energy equation: 

d{u6) d{v6) _ J. 


I I 

dt dx dy 


^ ^ ^ 

Pe V dy^ dz'^ 


(4.103) 


The initial and boundary conditions are: 


t = 0, 9 = 0 

(4.104) 

O 

11 

IcC: 

o' 

11 

(4.105) 

a; = 1, ^ = 0 

(4.106) 

= 0,| = -Nu 

(4.107) 

1! 

P- ' 

11 

o 

(4.108) 

„ 39 „ 

(4.109) 
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n 

^ = ^ = 0 (4.110) 

where, Nu is the local Nusselt nuinber. The rib is treated by the following con- 
dition. 


For X < li + and x > li, and y < drx 
dt Pe Oif Qy 2 q ^2 


For rr = or a: = -h and y < dri, kt 


dx 
39 

For X > or a; < + Iru and y = dri, ks tt 

dy 

4.3.2 Inverse Problem 


= k — 
solid ^ 


solid 


kf 


dy 


fluid 

fluid 


(4.111) 

(4.112) 

(4.113) 


The inverse problem is posed in the form of minimization of the following func- 
tional; tf 1 r 

J= [^ [ [ [(0^ - + w{9m - Ym)^] dxdzdt (4.114) 

Jo Jo Jo 

where and are the. computed and measured temperatures at the wall, Om 
and Ym are the computed and measured temperatures in. the interior and w is 
the weight assigned to the high quality data, 

4.3.3 Sensitivity Problem 

Sensitivity problem determines the the dependence of the temperature on the 
unknown local Nusselt numbers. To formulate the sensitivity problem, a per- 
turbation of Ad is given to 9 and 0 is replaced by 9 + A9 in Equation 4.103 to 
yield 


3(9 + A9) ^ d[u{9 + A9)] , d{vi9 + A9)]_ 

1- — ■+■ _ 


3t 

1 

Pe 


3x 


+ 


dy 


3Hd + A9) , 3^9 + A9) , 3^(9 + A9) 

1 — -|- ^ 


+ 


(4.115) 


3x^ ' dy^ ' 

Now, subtracting Equation 4.103 from Equation 4.115, we have the following 
sensitivity problem, 

N o/ A /i\ 1 va2l\Q\ Aff\ ^2/'A/9M 

(4.116) 


3{Ae) d{uA9) d{vA9) _ 

3x Pe 


52(A6') 3^{A9) d%A9) 

3x^ 3yY dz^ 
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In the initial and boundary conditions, Nu is replaced by Nu+ANu and 9 by 
9 + A^; thus we have: 


t = 0, 0 + A0 = O 

(4.117) 

(Z? = Oj 0 4“ A0 = 0 

(4.118) 

oy 

(4.119) 

!/ = o/(^ + ^^)= -(Nu + ANu) 

(4.120) 

y ~ d^ 9 A9 ~ 0 

(4.121) 

^ d{9 + A9) ^ 

.=0, ^=0 

(4.122) 


(4.123) 

Now, subtracting Equations (4.104-4.110) from Equations (4.117-4.123) respec- 
tively, we have the following boundary conditions for the sensitivity problem 

f = 0, A^ = 0 

(4.124) 

a; = 0, A9 = Q 

(4.125) 

o 

!l 

t-T 

11 

(4.126) 

„ d{A9) 

!/ = 0, g =-ANu 

(4.127) 

o 

11 

<1 

i 

II 

(4.128) 

-■T- 

(4.129) 


(4.130) 

The rib is treated by the following condition. 


For X < lx + Iri, and x > h, and y < dri 


d(A9) 1 as fd‘^{A9) ^ d^{A9) d^{A9)\ 

dt ~Feaf\ dx^ dy^ dz^ ) 

(4.131) 


. V 

For X = li ox X = lx Irii 2 .nd y < d^x, ks 
For X > or a; < lx + Itx-> V ~ ^ru 


= ki 


d{A9) 


solid 


dx 


g( Ag) 

solid ' 


fluid 

(4.132) 


fluid 

(4.133) 
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4.3.4 Adjoint Problem 


Let, A be the Lagrange multiplier that is associated with the constraint. We can 
write the functional 


J 



[{6^ - + wiBjn - Ymf']dxdydzdt + 

1 


A 


^ ^ ^ 
[Pe dx^ dy'^ 


d{u6) d{v6) dd 

dx dy dt 


dxdydzdt 

(4.134) 


In Equation 4.134, setting $ + A0 in place of 6 and J + A J in place of J, we 
have 

J + At7 = 


f [ [ [{9 + Ae-Y)^5{y)+w{9~b^0-Y)'^S(y-yi)]dxdydzdt + 
0 Jo Jo Jo 

tf f^d I^r [ d'^{e + AO) d‘^{9^A9) d'^{9 + A9) ' 

Pe dx"^ dy^ dz^ 


'0 



'0 


dxdydzdt 


r-tf 

>0 Jo 



A 


0 JO 


d{u{9 + A9)} d{v{9 + Ae)} d{9 + M] 

' dx dy dt 


dxdydzdt 

(4.135) 


where, 8{y - yo) i& the Dirac delta function centered around yo 

Subtracting Equation 4.134 from Equation 4.135 and neglecting the terms, 
involving squares of A9, we have 


AJ 


2 f ^ f [ f A9 [w{9 - y)%) + {9- Y)5{y - yi)] dydxdzdt + 
Jo Jo Jo Jo 

dxdydzdt — 



d‘^{A9) d^{A9) d\A9) 

dx^ + + dz^ 

d{u {A9)} ^ d{u(A^)} d{A9) 


dx 


dy 


dt 


dxdydzdt 


(4.136) 
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Using the boundary conditions from the sensitivity problem, we get 




, rs^) 

J dx^ 


dx 

d{A9) 


dx 


-1 a:=0 
■1 x^l 


x =0 


d{M) dX 
) dx dx 


X=l 

a :=0 


dx 


A9 


dx 


r a(A«)i 

1 1=0 1 


dx 


[ aa:J 


/'I «2 \ 

+ / 

J *=0 J 0 

/•I Q2\ 

+ / -;r-^{A9)dx 


x=l 


dx"^ 


( 4 . 137 ) 


Further 


ptf nl pd pr 

Jo Jo Jo Jo 


"*/ r d‘^(A9) 

A — dxdydzdt 


if fd r r 5(A^) 

A 


m r ' 

pf r<i r 

Jo Jo Jo 
pf p r 

Jo Jo Jo 


dx 


dx 


A9 


dydzdt 

J a:=0 

dydzdt + 


a;=l 


ftf pi pd pr q2\ 

AO-fr-r dxdydzdt 
ro dx^ 



( 4 . 138 ) 


Again 




fd^&o)^ ffax /■a“{A«)j 


dy"^ 


dy 

dm] 


dy 


j /==0 

y^d 


y-0 


f 


dm dx 

dy dy 


dy‘^ 
dy 


y—d 

y =0 


Jo 




Jy=0 


A9 


9y}y^0 

dyly^o 


I 


dy^ 

'^d^X 


+ / 


( 4 . 139 ) 
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Therefore 

ptf pi pd pr 

Jo 


,d^{A9) 

^ — dxdydzdt 

0 JO JO 



m r 

_ 

m r - 

_ . 


d{Ae) 


dy 


dy 


Ae 


y—d 

dxdzdt + 

y=0 

dxdzdt + 


J y=0 


if pi pd pr q 2 \ 

/ / / M^dxdydzdt 

0 Jo Jo Jo dy^ 


(4.140) 


Again 


- 


'0 


d‘^{Ae) 


J dz^ “ 

d{A9y 

Z^V 

dz 

2=0 

dX 

Z=T 

A0— 


dz_ 

2=0 


dX f d‘^{A6) 


^ f 

dzj 


dz"^ 


dz dz 


J 2=0 


dmd_x. 

dz dz 
d'^X 


=0 Jo dy^ 


(4.141) 


Therefore 

pf p p p ^ 

Jo Jo Jo Jo 


m d 


dz^ 


dxdydt + 


J z —0 


rtf p p r ^j^x, , , . 

/ / / / Ad-^dxdydzdt 

Jo Jo Jo Jo 


dz^ 


Also 


I 


^ ,d{u(A0)} 


A " = 

0 9x 




pl r)\ 

= [XuiAd)]lzl- 

= [AtiA^U - ^ u^{Ae)dx 


(4.142) 

"I X=^l • 

J x=0 

(4.143) 


Therefore, 


r^f p = Jp [XuAe]^^^dydzdt- 

Jo Jo 'do dO ■ 1 ' J ^ > 

jp u^{A9)dxdydzdt 


(4.144) 
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And 


I 




dy 


A 


/ 


d{^(Ag)} 

dy 


dy 


J \9y J dy 


1 y=:d 


dy 




y=0 


(4.145) 


Therefore 

f*tf nl nd nr 


ftf .i ra rr a{u(A0)} , , , r^f /•! fd pr Q. 

I I I I ^~^d.dydzdt^-J^ v-(Ae)d.dydzdt 

(4.146) 


And 


/" 


m. 


= [AA«]g' - r 

Jo 
pi} 

= A6^dt 


Ae^dt 

0 ot 

. M 


t=^tf 

t=0 


(4.147) 


Therefore 




tf pi pd pT 


ptj pL pa pr 

/ / / / A9-—dxdydzdt (4.148) 

Jo Jo Jo Jo at 

Now conabining Equations 4.136 to 4.148, we, have 


A J = 


pi fd pr^i /Q2X d'^x Q2x\ dx dX dX 
( A ^ + ttw + u — + + 


Pe Qy2 J 


dx 


dy dt 


Addxdydzdt + 


2 [^ [ f f [w{9 -Y)6{y) + i6- Y)5{y~ yi)]A6dxdydzdt 
Jo Jo Jo Jo 

f nA\M dyd.ii- f‘ r 

Jo Jo Vo L J a:=0 Jo Jo Jo Pe [ 

Pf p J ^A9 dxdydt— f f j [XuA9]^^^dydzdt 
Jo Jo /o Pe [ay 1^=0 •/o Jo Jo 

m [XA9]i^i dxdydz 

' 


dydzdt 


(4.149) 
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are 


For A J ->■ 0, the following adjoint equation and final boundary conditions 
deduced 


t — tf, A = 0 

x = 0, A = 0 

ic = 1, -T — h uAPe = 0 
ox 

o "V 

y = 0, = -2u;Pe(0i„ - Fu,) 

y = d, A = 0 

" = “• a; = “ 

n 


The rib is treated by the following condition. 


For x<li + Iri, and x > h, and y < dri 
m~ Feaf\dx^ dy^ dzy 


For X ~ 1% ov X = h Iri) and y < dri, 


Since 


1 r 


dX 


= kf 

dX 

jU 

ox 

solid 

dx 



= kf 

dX 

ks ^ 
dy 

solid 

dy 

-ANu)]y^Q 

dxdzdt 


fluid 

fluid 


(4.151) 

(4.152) 

(4.153) 

(4.154) 

(4.155) 

(4.156) 

(4.157) 


(4.158) 

(4.159) 

(4.160) 

(4.161) 


(4.162) 


and the definition of a gradient 

r J'{-mvi)dxdzdt 

Jo Jo Jo 

From Equations 4.161 and 4.162, we have the following expression for the 
gradient of the functional. 

A J _ A(t, xA^) (4.163) 


J' = lim , T 

aNu-+o ANu 


Pe 
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4.3.5 Determination of Search Step Size. 

If the value of the functional is after k-th. iteration and it becomes after 
temperatures are updated through a step size of yd*, the condition for the optimum 
step size will be given by 


djB^ 


(4.164) 


Recalling the expression of functional 

f*tf pi pd pr 


pif pi pa pr 

~ Jo Jo Jo Jo ~ ~ (4.165) 


Therefore 


mw^J*(-AAfn)* 


minp r f r r{[0{(-AiVn)* - ^*P*} - Yf5{y - y{) 
Jq «/o ^0 

+n;[6'{(- AiV«)*= - /3*P*} - Yf5{y)}dxdydzdt 

(4.166) 


linearizing by Taylor’s series expansion 


j*+i[(_AArn)*+^] = 

r f f f p[(-AArti)* - yd*A0(P*)] - F]^ 5{y - yi)dxdydzdt 

t/o t/n Jq Jq 


fO t/o ^0 JQ 

rtf pi pd, pr 


+W f f f [0[(-AAru)* - ^*A0(P*)] - yf 5{y)dxdydzdt 
Jq Jq Jq Jq 

f f { [{d~-Y)-P’'A9y{5{y-yi)+w5{y)]dxdydzdt 
Jq Jq Jq Jq 


(4.167) 


Therefore 

f f f A0{9-Y-^'°Ae){5{y~yi)+w6{y)}dxdydzdt (4.168) 

Jo Jo Jo Jo 

Using Equations 4.164 and 4.167, we have the following equation for the op- 
timum step size 


k _ So fJ fJ Jo ^ ~ + wS(y)} dxdydzdt 

^ ^ /o^ Sj So [^(y ~ 2/i) + w5 (y)] dxdydzdt 


(4.169) 
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4.3.6 Complete Inverse Algorithm. 

The complete inverse algorithm for calculation the Nusselt number distribution 
over flat surface is summarized below; 

1. Solve continuity and momentum equations to get the velocity field. 

2. Smooth imperfect plate temperature data using Gram orthogonal polyno- ' 
mial filter. 

3. Solve energy equation using the following initial and boundary conditions: 


t = 0, 0 = 0 

(4.170) 

a: = 0, 0 = 0 

(4.171) 

dO 

(4.172) 

y = 0, 9 dw 

(4.173) 

II 

11 

o 

(4.174) 

hi 

II 

o 

Cd| 

hi 1 

II 

o 

(4.175) 

n 

(4.176) 

The rib is treated by the following condition. 


x<li + Iri, and X > k, and y < dri, u = 0, v = 0 

(4.177) 

G6t the vsblues of loc3.1 Nusselt numbers. Tnke an initial guess of Nusselt 
number to be same as that found at the end point of the domain. Initialize 

J = 0 


Solve energy equation using the following initial and boundary conditions 

(with Nusselt number from step 4): 

O 

il 

C> 

11 

(4.178) 

a? == 0, 0 = 0 

(4.179) 


(4.180) 
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de 

y = 0, — = -Nu 
oy 

y = d, 6 = 0 

. = 0 , 1 = 0 


89 

dz 


z = r, — = 0 


The rib is treated by the following condition. 

X <li + Iri, and X > h, and y < dri, u = 0, v = 0 


(4.181) 

(4.182) 

(4.183) 

(4.184) 

(4.185) 


Find out the temperatures at the plate as well as the plane of measurement 
of redundant data. 


6. Compute the functional, 

J = / ^ f f {{Oyj - + w{9ra - Ymf}dxdydzdi 

Jo Jo Jo Jo 

7. Check the stopping criterion, 

where, e is a predetermined real number, very close to zero. 

If it is satisfied, stop iterations. Otherwise, 

8. Solve the adjoint problem 


(4.186) 


(4.187) 


Pe 


X = 1 , 

n 

= 0, ^ 

dy 


- + - r)% 

dy at 

- yi) = 0 

(4.188) 

ztfyX = 0 

(4.189) 

= 0, A-0 

(4.190) 

|^4.^iAPe=:0 

ax 

(4.191) 

= -2ioPe(0u; - Yyi) 

(4.192) 

= d, A = 0 

(4.193) 

n n 

(4.194) 
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The rib is treated by the following condition, 

X <h + Iri, and X > k, and y <dri, u~ 0 , -u = 0 

9 . Compute the gradient 


J' = 


A(t,x,0, z) 

pi 


10. Compute the descent direction 

pk ^ ^kpk-i ^ j,k 

where 

for first iteration 7* = 0 
otherwise 7'' = 




(4.195) 

(4.196) 

(4.197) 

(4.198) 

(4.199) 

(4.200) 


11. Solve the sensitivity problem 


d{Ae) , d{uA9) , divAO) _ 1 


d'^jAe) d‘^{Ae) d^jAe) 

dy Pe [ dx^ dy"^ dz^ 


dt ' dx 

with the following initial and boundary conditions; 


(4.201) 


X 


y - 0, 


It 

0 

t> 

II 

0 

(4.202) 

0 

If 

<I 

o' 

11 

(4.203) 

= 1 

’ dx 

(4.204) 

= -ANu = 
dy 

(4.205) 

y = d,Ae = 0 

(4.206) 


(4.207) 


(4.208) 
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12. Compute the step size 

fo^ fo So ~ - Vj) + Hv) dxdzdt 

p = ^ — J (4 209) 

/o' /o' /o' [ES «il{v - Bi) + -S(!/)l dxdzdt 


where mp indicates number of measurement planes and Wj is the weight 
assigned to the data from the j-th plane of measurement. 


13. Compute new local Nusselt number 

(_Nu)^+^ = (-Nuf - 


(4.210) 


14. Go to 5 and continue till convergence. 


4.4 Transient Plate Temperature Problem 

The inverse formulation for an transient surface temperature is studied here. Both 
the velocity and temperature fields are two dimensional. No additional interior 
measurement is considered. 

4.4.1 Assumptions 

1. Temperature field is two dimensional, i.e., plate temperature is varying with 
length only. 

2. Plate temperature is unsteady. 

3. The local Nusselt number distribution is not varying with time. 

The local Nusselt Number is defined as 

Nu(x) = — ■ (4.211) 

4.4.2 Direct Problem 

The direct problem solves the following energy equation: 

de d{u6) ^ d{v6) _ _r 
dt^ dx dy Pe dy^ J 


(4.212) 
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The initial and boundary conditions are: 


i = 0, ^ = 0 

a; = 0, 1^ = 0 

x = i, e = Q 

2/ = 0, = -0Nu 

ay 

y = d, 9 = 0 


(4.213) 

(4.214) 

(4.215) 

(4.216) 

(4.217) 


where, Nu represents local Nusselt number. The rib is treated by the following 
condition. 


For X < li + Iri, and x > li, and y < dri 
dd 1 o, /dH d^\ 


dt Pea/ \dx^ dy'^ 


(4.218) 


dd dO 

For X h ov X h + lri> and y < dri, k, ~ —^f'^ (4,219) 

solid fluid 

do dO 

For X > h ox X Iri, and y = d^i, k, — = kf — (4.220) 

% I solid ^1/ fluid 

4.4.3 Inverse Problem 

The inverse problem is posed in the form of minimization of the following func- 


tional: 


J= f f [(9- Y)^]6{y)dxdydt (4.221) 

Jo Jo Jo 

where ^ and Y are the computed and measured temperatures and 5 is the Dirac 


delta function. 


4.4.4 Sensitivity Problem 

Sensitivity problem determines the dependence of the temperature on the un- 
known local Nusselt numbers. To formulate the sensitivity problem, a pertur- 
bation of A9 is given to ^ and 9 is replaced by 9 + A9 in Equation 4.212 to 


yield 
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Now, subtracting equation(4.212) from equation(4.222), we have the following 
sensitivity problem, 


d{Ld) djuLO) d{vl\B) 1 

dt dx dy Pe 


d^{A9) a^(Ae) 

dx^ dy'^ J 


(4.223) 


In the initial and boundary conditions, Nu is replaced by Nu+ANu and 9 by 
9 + A9, we have: 

t = 0, 0 + A0 = 0 (4.224) 

x = 0, ^ + A^ = 0 ■ (4.225) 

(4.226) 
oy 

d{9 + A9) 


dy 


y = 0, 


= -(Nu + ANu)(0 + A0) 


(4.227) , 

(4.228) 


dy 

y — dy 9 A9 — . 0 

Now, subtracting Equations (4.213-4.217) from Equations (4.224-4.228) respec- 
tively, we get the following boundary conditions for the sensitivity problem: 


t = 0, A9 = 0 

x = 0, A0 = O 


y = 0, 


X 

d{Ae) 


= -9ANu - NuA^ 


dy 

y = dy A9 = 0 


(4.229) 

(4.230) 

(4.231) 

(4.232) 

(4.233) 


The rib is treated by the following condition. 


For X < li + Iri, and x > li, and y < dri 

a(A0)_J_^' 


a^jAg) a^(Ag) 
arc^ at/2 


d{A9) 


ai Peo!/ 

For a: = or a: = + ^rii and y < dri, K 

d{Aff) 


d{A9) 


For x>li ov X < and y dri, ks 


dy 


solid 

isoud 


(4.234) 


flmd 

(4.235) 


fluid ^ 
(4.236 
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4.4.5 Adjoint Problem 


Let, A be the Lagrange multiplier that is associated with the constraint. We can 
write the functional. 

rtf rl rd 


j 


-L 



{9 — Yy5{y)dxdydt + 


0 Jo Jo 

if pi pd 



X 


'0 Jo Jo 


^ d{u9) d{v9) 89 


[Pe \dx‘^ dy'^ ) dx dy dt 


dxdydt 

(4.237) 


In Equation 4.237, setting 9 + A9 in place of 9 and J + A J in place of J, we 
have 

ftf pi pd 


J + AJ = [ [ [ {9 + A9 -YY6{y)dxdydt + 

Jo Jo Jo 

A 

Pe 

nl pd 


d^{ 9 + A9) d^{9 + A9) 

dx"^ dy^ 


dxdydt ■ 


m o 


A 


d{u{9 + A9)} ^ d{v{9 + A0)} ^ d{9 + A9) 


dx 


dy 


+ 


dt 


dxdydt 

(4.238) 


Subtracting Equation 4.237 from Equation 4.238 and neglecting the terms 
involving squares of A9, we have 

f'tf pi pd 


AJ = 2 [ ^ [ [ A9{9 ~-Y)5{y)dxdydt + 

Jo Jo Jo 

Jo Jo Jo P® 


dxdydt 


m d 

_ 


A 


d^A9) d^{A9) 

dx"^ dy"^ J 

9HA0)} , a{u(Ag)} , d{A9) 
dx dy dt 


dxdydt (4.239) 


Using the boundary conditions from the sensitivity problem, we get 

r fd\A9), frdXfdHm.\. 

' / 


= 

^ dx^ 


J dx"^ 

d{A9y 

x=l 

dx 

x^O 

d{A9y 

x=i 

dx _ 

a;=0 

r a(Ae)i 

A o 


X—l 

a;=0 


dx dx 

rlQ2x 


dx 


dx 


aj—O 


z;=:4 Pi 

- A«|^l + [‘ A(A«)<fa 
L dx j^^i Jo dx^ 


(4.240) 
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Further 


r = - f 

Jo 



if pd nr 


fQ Jo Jo 


dx^ 



A 


0 JO JO L 

tf pd 


a(Ag) 

dx 


dydt 


n . 

i 


dx 


M 


J x=0 

dydt + 


-i X=1 


if rl nd q2^ 

d.d~dxdydt 



(4.241) 


0 JO JO 


Again 




d'^jAe) 
dy'^ ^ 

d{Ae)'^^=‘^ 


dy J 

d{Ae) 

9y Jj,=o 


-/(I/ 

I. 


d^AO ) 
dy^ 


dy dy 


-i j/=d 
2/=0 


djAd) dx 
y=o Jo dy 9y 


dy 


iy=d 


n3/=d 

dy 


A^l + 


A6 

A6 


dy. 

dX 


y-d 

,.aL 


(4.242) 


Therefore 


if f‘ ^ d\A$) 
dy^ 


fif n r 

Jo Jo Jo 


X- 


-dxdydt 


■ 

_ 

n i 

_ 


d{A9) 
^ dy 


dy 


A9 


y=d 

dxdt + 

y=0 

dxdt + 


Jy=0 

m d f)2\ 


(4;243) 


Also 


L 


' , 5{u(A0)} 


A 


dx 


dx = 


A 


/ 


g{u(A^)} 

dx 


dx — 


dx f d{u{Ae)y 


dx 


I 


dx 


dx 


X=J 


x=0 


(4.244) 


Therefore 


if 


r^dx 

= [A«(Ae)i‘:|, - 1 

ydX 

= [A«A4=, - 

f\ f X?Jyj^i^^dxdydt ^ [^ [ [XuAei^^dydt- 
"" ' %f h dx 

/ / u^{Ae)dxdydt {4.24b) 

Jo Jo Jo 
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And 




dy 


A 


I 




dy 


/(!/ 


dy 


dy 


yz=d 


J 2/=0 


rd 

= [Hmizt - 1 

I'd f)\ 

= [A»(A9)],,, - v^mdy 


(4.246) 


Therefore 


r [' = - r f rv^^md.dydt 

Jo Jo Jo % JO Jo Jo dy 

■ [^ [ [Xv{/:\9)]y^^dxdt (4.247) 

Jo Jo 


And 


L 


. d{Ad) 
^ dt 


dt 


+ 




dx 


dt 

= lAAC5 '-/‘'Ae| 
= lAA0L,-f 


dt 


AOdt 
dt 
dt 


t=0 


(4.248) 


Further 


l\^^dxdydt = 1 1 

(4.249) 


Now combining Equations 4.239-4.249, we have 


AJ 


m d 

^ ■ I d 

2 f f {d-Y)6{y)A6dxdydt- li [XAe]t^t^dxdt 

Jo Jo Jo . U r r^\l 

r- fAhm] dyM- r 


L(^ + ^A+ui^ + v^ + 

Pe \9x^ V 


8A , ex ^ ax 
dx dy at\ 


Addxdydt + 


J aj=0 


r*/ 1 

Pe 


A9 


dx 


dydt — 


?’ f [a»a«!,=, iydt - jj I [A»(A«)1,=. 

Jo Jo J ^ r r^\ 

f’l 


J a:=i 


-(- 


■*/ 7'_1_ 

Pe 


.a(Ag) 
dy J 


dxdt + 


/*/ 

7o io ^ 


% 


dxdt (4.250) 


J2/=0 
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For A J — >■ 0, the following adjoint equation and final boundary conditions are 
deduced: 


Fe\dx^^ dy^)'^ dx ^~Tr 
t = tf, A = 0 

a; = 0, A = 0 
dX 

X = 1, — + uAPe = 0 
ox 

o 

y = 0 , ~ + ANu + 2 Pe (0 - F) = 0 
y = d, A = 0 

The rib is treated by the following condition. 

For X < li + Iri, and x > li, and y < d^i 

^ _ __L£i 

dt Pea/ dy^ J 


Since 


For a; = /i or a: = /i + 4i) and y < dru K 
For x>li ox X <h + ^ri) and y = d^, ks 

1 


dx 

dy 


solid 


kf 


dx 

= k — 
solid ^ 


fluid 

fluid 


(4.251) 

(4.252) , 

(4.253) 

(4.254) 

(4.255) 

(4.256) 


(4.257) 

(4.258) 

(4.259) 

(4.260) 


[Aj^(ANu)]j,^o dxdt 

and the definition of a gradient 

A J = / ^ / J'(ANu)da:di 

Jo Jo 

From Equations 4.260 and 4.261, we have the following expression for the 
gradient of the functional: 


(4.261) 


AJ 

/(A = lim -T-ry- 
aNu-40^Nu 


ftf 

Jo 


xe 

Pe 


di 


(4.262) 


2/=0 
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4.4.6 Determination of Search Step Size 

If the value of the functional is after A;-th iteration and it becomes after 
temperatures are updated through a step size of the condition for the optimum 


step size will be given by 




= 0 


(4.263) 


Recalling the expression of functional 

rtf pi rd 


~ lo lo lo (4.264) 


Therefore 


minpj'^i-ANu)'^ = min^ f f f [^{(-AiVu)^ - - Ff 5(y)da:dydt 

Jo Jo Jo 

(4.265) 


linearizing by Taylor series expansion 

rtf rl I'd 

>0 Jo Jo 


jfc+i[(_AiVu)''+i] = [^ f [ [^[(-AiVu)^ - ;0*A^(P*=)] - yY 5iy)dxdydt 

Jo Jo Jo 

(4.266) 

Therefore 

Qjk+l ptf pi pd 




2 f f Ae{e - Y - P^Ae)5{y)dxdydt (4.267) 

Jo Jq Jo 


Using Equations 4.263 and 4.267, we have the following equation for the op- 
timum step size 

j:‘j!.fAs{s-my)d^dvdt 

lo' 


(4.268) 


4.4.7 Stopping Criterion 

The additional measurement approach [Section 2.5.2] is considered as the stopping 
criterion. An additional functional Ji is constructed with extra temperature data 
that is known to have small errors. The iteration stops when Ji starts to increase. 

4.4.8 Complete Algorithm 

1. Solve continuity and momentum equations to get the velocity field. 

2. Take an initial guess of Nusselt number. Initialize J ^O 
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3. Solve energy equation using the following initial and boundary conditions: 


t = o, e = o 
x = 0, 0 = 0 


X 


90 „ 

1 , — = 0 

ox 


v = (i. = 

y = d, 0 = 0 

The rib is treated by the following condition. 

X Iri, and X > li, and y < « = 0, v = 0 

Find out the temperatures at the plate. 


(4.269) 

(4.270) 

(4.271) 

(4.272) 

(4.273) 

(4.274) 


4. Compute functional, 

J = f f [{0- F)^] 6(y)dxdydt 
Jo Jo Jo 

5. Check stopping criterion 

where, e is a predetermined real number, very close to zero. 
If it is satisfied, stop iterations. Otherwise, 

6. Solve adjoint problem 

1 Z'gA d{uX) ^ djvX) dX 

Pe \dx^ '^ dy^ J 9x dy dt 

t — tftX = 0 

X = 0, A = 0 

X = 1, uXPe = 0 
ox 

y — Q — + ANu + 2Pe(^ ~ 
dy 

y = d, X = 0 

The rib is treated by the following condition. 

X < li + Iri, and X > h, and y < dri, u = 0, v = 0 


(4.275) 

(4.276) 


(4.277) 

(4.278) 

(4.279) 

(4.280) 

(4.281) 

(4.282) 

(4.283) 
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7. Compute the, gradient 




\xe 


Pe 


dt 


Jy=0 


8 . Compute the descent direction 


where 


pk ^ ^kpk-i ^ j,k 

for first iteration 7 *^ = 0 

(jfc_ jfc-i|jfc) 


otherwise 7 * = 




9 . Solve the sensitivity problem 


djAO) d{uA9) d{vA9) ^ 1 

dt dx dy Pe [ dx'^ 


d\A9) d'^{A9) 

" 1 “ 




with the following initial and boundary conditions: 

i = 0, = 0 


a; = 0, A9 = 0 
d{A9) 


X = 1, 


dx 


Q 


y = Q, 


10. Compute the step size 


d{A9) 


dy 


-f A^Nu = -ANu0 = P'‘e 

= d, A9 = 0 




n f:A9{e -Y)5{y)dxdydt 
j:PiJ^iA9f5(y)dxdydt~^ 

11 . Compute the new local Nusselt number 

(-Nu)*+^ = (-Nu)* - 


(4.284) 


(4.285) 

(4.286) 

(4.287) 


(4.288) 

(4.289) 

(4.290) 

(4.291) 

(4.292) 

(4.293) 


(4.294) 


(4.295) 


12. Go to 3 and continue till convergence. 



Chapter 5 

Apparatus and Instrumentation 


The present study is concerned with the measurement of heat transfer coefficient 
and the Nusselt number a heated surface of a rectangular duct exposed to a flow 
field. Figure 5.1 shows the sketch of the experimental setup used in this work. It is 
identical to the one used by Singh [2002]. The experimental facility comprises of a 
flow circuit (wind tunnel), the heating section, traverse mechanism and an image 
processing system. Hot wire anemometry (HWA) and resistance thermometry 
(RTD) have been utilized for velocity and temperature measurements adjacent 
to the surface carrying the rib. Local temperature of the plate has been measured 
using K-type thermocouples connected to a National Instrument data acquisition 
card (NI 4351). Wall temperature distributions has also been recorded with liquid 
crystal sheets, exposed to two 150W tungsten-halogen lamps. The lamps provides 
an excellent color rendering even in the long run, a reasonable high efficiency, high 
luminance and compact size, as described by de Boer and Fischer [1981]. 

5.1 Wind Tunnel 

The wind tunnel is of the open loop type which operates with the fan in the 
suction mode. The ambient air flow is sucked from the temperature controlled 
room into the test section through a flow straightener and five screens in the 
settling chamber and a 3 ;1 contraction cone. Air flows through the bell-shaped 
contraction cone, the heated section, the unheated length of the duct and is then 
exhausted by a centrifugal fan that is run by a 3-phase motor.The speed of the 
blower is controlled by a speed controller (Victor G 1000) supplied by Kirloskar 
Electric Co , India The settling chamber is 1950 mm long and has a rectangular 



5.2 Heating Section 


92 



Figure 5.1: Schematic drawing of flow system, coordinate system and instrumen- 
tation. 

cross section of 500 mm by 1000 mm. The test channel is 3300 mm long with an 
aspect ratio of 1.8:1, (160 mm x 298 mm in the vertical plane). The test channel 
is made of perspex sheet of 12 mm thickness. Velocities in the range 0.3-4 m/s 
are presently realizable in the tunnel. The turbulence intensity in the incoming 
flow was found to be less than 0.5%. 

5.2 Heating Section 

The heating section is an assembly of an aluminum plate, Bakelite sheets, and a 
stainless steel foil. 

The heating section is situated just at the entrance of test duct. The test plate 
is 730 mm long and is located flush with the bottom wall. It serves as the heat 
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transfer surface, both, as smooth and one roughened by a rib, while the top and 
two vertical walls are smooth and thermally insulated. Downstream of the test 
duct is an unheated section made of perspex with the same cross section. The 
heat transfer surface is constructed using a single aluminum plate (680 mmx 298 
mm X 3 mm). There are six stainless stainless steel foils of dimension 690 mmx 47 
mmx 0.045 mm that are cemented on to the 25 mm thick Bakelite sheet. These 
foils are electrically connected in series and finally connected to a DC source for 
power supply. To prevent the aluminum plate from electrical contact the upper 
surface of the stainless steel foil is coated with varnish. A 2 mm mica sheet is 
placed between the Bakelite and aluminum plate for thermal insulation. The 
surface of the aluminum plate is highly polished to minimize thermal emission, 
and hence the radiative losses. In addition, to minimize the conductive heat 
losses, the lower surface of the Bakelite board is insulated using a 13 mm thick 
of Bakelite with a 2 mm thick air gap between them. 

The heat transfer surface is also instrumented with thirteen calibrated ther- 
mocouples of the chromel-alumel type. The thermocouples are located inside 
the heated surface with omega-bond 101 epoxy. Twelve thermocouples are dis- 
tributed along the center-line of the heated plate for wall temperature measure- 
ment. In order to check the span-wise temperature uniformity, surface tempera- 
ture along the lateral of the heated plate is also instrumented with one thermo- 
couple. To maintain a smooth flow surface, the thermocouple beads axe attached 
to the underside of the aluminum plate fed through the opening made in the two 
Bakelite sheets and mica. To measure the conduction losses on the underside 
there are two thermocouple are mounted on the two sides of the upper Bakelite 
sheet. There are two more thermocouples mounted at the two sides of the lower 
Bakelite sheet. The four thermocouples fall on one line. This arrangement is ■ 
repeated at a second location in the downstream. 

5.3 Hot-Wire Anemometry 

The hot-wire anemometer is an example of an intrusive technique whose appli- 
cation for measuring turbulent flow has far outstripped other mstruments The 
hot-wire is basically a thermal transducer and responds primarily to the velocrty 
magnitude and is based on the principle of heat transfer. This sensor is small rn 
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size and coupled with a feedback circuit, it has a very high frequency response. 
The DANTEC model 56 C 17 CTA has been used in the present study. The main 
unit 56 C 17 CTA delivers the servo-voltage as the output of the instrument. This 
voltage is a measure of the fluid velocity. The feedback circuit of the CTA plays 
an important role in improving the frequency response of the hot-wire from about 
lOOHz to lOKHz. 

The 56 C 01 CTA contains a function switch with three modes for operation, 
namely TEMP, STD.BY and FLOW. In TEMP position the resistance of the 
connected probe can be measured in terms of a current supplied to it. In STD.BY 
position no current flows through the bridge. In FLOW setting the CTA starts 
operating with the function of the servo amplifier. A setting named BRIDGE ADJ 
enables the adjustment of bridge balance for measurement of probe resistance and 
setting of the desired overheat resistance. In addition, BRIDGE ADJ has a switch 
pair for coarser adjustment of overheat resistance and a screw for fine adjustment. 
Resistance settings ranging from 0-30 in steps of O.OOIQ are possible. This 
adjustment is crucial for adjusting the overheat resistance for the calibration 
procedure. The CTA in TEMP mode produces a voltage proportional to the 
resistance of wire. The Calibration of this voltage against a known temperature 
can be used to measure the instantaneous temperature of the probe. The mean 
value unit 56 N 22 is a 5.5 digit display voltmeter. The primary purpose of this 
module is to measure the DC component of the output signal from 56 C 17 CTA. 
This module has 100 /ivolt resolution, 1-1000 second integration time and 14 
switch selectable inputs, 

The accuracy of hot-wire measurements is affected by the accuracy of the 
calibration procedure that of the procedure that is used to solve the nonhnear 
simultaneous equation. Calibration of a hot-wire probe involves two major steps, 
namely data generation and curve fitting. Calibration data is generated by mea- 
suring the output of the anemometer when the probe is subjected to flow with 
known velocity. In the present study, calibration has been performed in the test 
cell itself due to its high quality uniform inflow. 
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5.4 Resistance Thermometry 


DANTEC Platinum coated tungsten wire ( 0:20 = 0.36/°C) probe (Pll) has been 
utilized as a temperature measuring device. The resistance of the wire of the 
probe at 20'’C has been determined by a 5.5 digit display HP 3457A multimeter 
against the precalibrated thermister TBX-68T mounted on the terminal block of 
the National Instrument’s data acquisition card (NI 4351). The HP 3457A has 
a reading and program storage module and is available of making unsteady mea- 
surements. The HP 3457A’s math operations manipulate or modify a measured 
reading before it is displayed; in addition the STAT operation evaluates mean 
and standard deviation. HP 3457A is programmed for converting the resistance 
in terms of the operating wire temperature. The temperature of the fluid flow is 
measured by using the relation: 

RiTsensor) — -^20 P 0:20-^20 (Tsensor ~ T 20 ) (5-1) 

The RTD probe is attached with a computer controlled traverse mechanism 
for accurate positioning. 


5.5 Liquid Crystal Thermography (LCT) 

The liquid crystal technology is one of the important inventions of the twentieth 
century that has found many applications as a visualization tool today. Liq- 
uid crystals change their color on applications of external stimulus, for example 
temperature and shear stress distributions and thus act as a measure of their 

change. 

Liquid crystal is a unique substance, which exists between the solid and the 
isotropic-liquid phase of some organic compounds. It scatters incident light very 
selectively. The lowest temperature where liquid crystals scatters visible light 
is called the event temperature. At a temperature below the liquid crystal's 
event temperature, the liquid crystal wiU be in the solid state and will appear 
transparent. The clearing point temperature is the temperature at whi^ t^e 
liquid crystal ceases to reflect visible light. When the temperature ^ceeds tte 
clearing temperature point, the liquid crystal will enter the pure liquid state 
and appear transparent. The color output in liquid crystals is reversible and 

reproducible. 
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5.5.1 Forms of Liquid Crystals 

Liquid crystals are commonly categorized based on their molecular morphology 
and optical properties, e.g. smectic, nematic, and cholesteric. Cholestric LCs, so 
named because the first of this type were derived from cholesterol, are capable 
of displaying a range of visible colors and are thus classified as thermochromic. 
However non-sterol based thermochromic LCs have also been developed, and are 
termed chiral-nematic, in reference to their molecular arrangement. 

The nematic structure comes about when long range inter-molecular forces 
cause molecules to align relative to one another in planes. These planes are de- 
termined by the local solid boundaries and the direction of alignment of molecules 
within a plane is called director. Direction of this director changes on the ap- 
plication of mechanical or electrical deformation. The cholesteric materials form 
a similar structure. However, due to the shape of the molecules, the director of 
adjacent planes are rotated through certain angle and thus forming a helix of a 
certain pitch. The pitch of the helix changes depending upon the temperature of 
the liquid crystal. For the liquid crystal cholesteric materials considered for tem- 
perature sensing purpose, It is necessary that the pitch of the helix represented 
by the director is of the order of the wavelength of the visible light. 

5.5.2 TLC Surface Coating 

The most common means of applying thermochromic liquid crystals to a surface 
is by using micro-encapsulated form of the material. The tiny spherical capsules 
(20/xm in diameter) are glued to the test surface using an adhesive. Care is taken 
to ensure that the surface is fully dry before the mixture of crystal and binder 
mixture is applied. A thin (20/i mm thick) coating is best developed by spraying 
many thick coats, with interim drying using an air heater. It is also worth noting 
that thin plastic sheets are available fi:om suppliers and these often display a very 
clear color display. TLGs are available in different bandwidths, which defines the 
temperature range in which it actively reflects visible light. In the present work 
the commercially available thin plastic sheet of TLC (Hallcrest, Inc.;R35C5W) is 
pasted on the heating surface. This means the activation of the red color of this 
particular TLC begins at 35‘’G and the bandwidth is 5 degrees. The advantage 
of a precoated sheet are convenience and quality of color but these must be set 
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against possible difficulties in wrapping the sheet over complex surfaces with 
compound curvature. This is not the case for the present work. 

5.5.3 Response Time 

During the transient heat transfer experiment, the accuracy of the measured 
temperature depends on the rate at which the liquid crystal’s optical properties 
respond to the change in the surface temperature. Ireland and Jones [1987] 
measured the response of nematic type LC material and they were able to quantify 
the lag in the optical signals at rate of increase of surface temperature as high as 
2000 °C/sec. Experiments showed that the delay between the time at which the 
surface reached a steady state color display temperature and the occurrence of the 
color display was no more than a few milliseconds (typically 3 to 10 milliseconds). 
In fact, the time taken for any point within the liquid crystal layer to achieve 
the heated surface temperature is a function of the distance from the surface, the 
layer thickness and the film diffusivity. 

5.5.4 Range of Application and Limitation 

Although the TLCs offer an excellent way of measurement of full surface temper- 
atures there are always some constraints on their use. 

Temperature Range 

Encapsulated liquid crystal materials are available in the temperature range -30 
to 115 °C. The temperature range over which a material is optically active (the 
color play) can also be selected. The narrowest color range is 1.0 -C and the 
widest is about 20 "C. Wide-hand crystals typicaUy have bandwidths between 
5“C and 20°C. They are useful when an object has large temperature variations. 
Practical applications of wide-hand crystals include investigation of the surface 
temperature distribution on gas turbine blades or characterizing the temperature 
distribution on electronic components. 

Pressure Range 

Encapsulated TLCs are insensitive to pressure. They are tested up to 133 bar 

without any appreciable change in its behavior. 
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Speed of response 

As discussed earlier, the optical response time of TLCs are in the range of few 
milliseconds e.g. for chiral-nematic it is around 3 msec. 

Chemical Contamination 

TLCs are complex organic compounds that are particularly susceptible to attack 
from organic solvents. In addition, UV (ultraviolet) light can destroy the color 
display, so care must be taken to minimize exposure to UV light. 

Maximum Heat Flux 

As a rule-of -thumb, a value for maximum heat flux of 2 x 10^ KW/m^ can be 
taken. The conductivity of the liquid crystal layer is typically 0.2 W/mK and 
hence the above heat flux gives a 2^0 temperature drop across a 20 pm layer. 

5.5.5 Image Processing 

In a transient experiment, any site on the TLC surface will change color as the 
surface passes through optically active temperature range. This permits either 
color processing or intensity based (black and white) processing of the image of 
the surface optical response. 

The standard color television practice of representing color by R, G and B 
tristimulus signals, red, green and blue, is used in color analysis. The tristimulus 
signals can be thought of as outputs from three camera detectors each with its op- 
tically sensitive range centered on a different wavelength in the visible spectrum. 
The intensity signals is made up of the sum of the R, G, and B signals. The R, G, 
and B signals are changed to H, S, and I (namely hue, saturation and intensity) 
form, which is less sensitive to the variations in strength of illumination. The 
derived signals hue and saturation are both independent of intensity and are pro- 
duced using look-up tables in commercially available color frame grabbers. Hue 
depends on the ratio of RGB signals and its use has found favor in TLC work 
since its value normally increases monotonically with temperature. Hue physi- 
cally represents the dominant wavelength of the light being displayed by the LGs. 
It is determined by establishing the angle between the orthogonal red, green, and 
blue components in RGB space. 
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Various techniques have been employed to calculate hue from the RGB com- 
ponents of the digitized images. The hue is calculated as (Ireland and Jones, 
1987): 

2 Q ^ 

“ {6[i? - 7)2 + (G - 7)2 + (B- /)2]}o.5 

where I is defined as: 


R G B 

3 


(5.3) 


The value of hue varies from 0 to 180°. The highest value of hue corresponds to 
perfect blue while the lowest is red. 


5.5.6 True Color Image Processing System 

The image acquisition and processing system used in the present investigation • 
consists of a high resolution 768 (horizontal) x 574 (vertical) pixels CCD video 
camera (SONY XC-003P) with 16 mm focal length lens (VCL-16WM), a 24-bit 
color frame grabber board (Imaging Technology) and a high speed PC.These sys- 
tems store the appropriate intensities of red, green, and blue needed to produce a 
corresponding matched color response at each point in an image. The frame grab- 
ber can be programmed for color analysis using base level C-programs. Several 
software programs are commercially available to analyze images. The software is 
capable of acquiring, freezing, sequential grabbing of images with specific func- 
tions such as determination of red, green, blue values and also the corresponding 
hue distribution. The current 24-bit true color image processing board from 
Imaging Technology, USA is programmed for the aforementioned color analysis 

using base level C-programs. 

5.6 Uncertainty and Measurement Errors 

Possible sources of errors in the experimental data ate; (a) fluctuations m the 
supply voltage of the D.C. power source, (b) the positional accuracy in locatmg 
the probe, (c) errors in calibration data, (d) inadequate compensation tor room 
temperature, (e) inadequate signal length and sampling rate, (t) dnft in elec- 
tronics and (g) non-uniformity in Ulumination on the liquid crystals. Most of 
the experiments have been conducted several times and the repeatabihty of the 
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results presented have been confirmed. The local Nusselt number and average 
Nusselt number were found to be in good agreement with published correlations 
and energy balance check [Tariq et al, in press]. 


Chapter 6 

Data Reduction 


Data reduction for velocity, temperature and the local Nusselt number is discussed 
in the present chapter 

6.1 Velocity Data Reduction 

One dimensional velocity data was collected using a hot-wire anemometer. A 
single wire probe was used for this purpose. The data obtained was not in terms 
of velocity but it was in terms of voltage which is required in excess to maintained 
the wire at its original temperature. There is a unique relationship between this 
voltage and velocity, and the relationship must be obtained by proper calibration. 

6.1.1 Hot-wire Calibration 

Calibration of a hot-wire probe involves two steps, namely data-generation and 
curve-fitting. 

Data-Generation 

Calibration data is generated by measuring the output of the wire when it is sub- 
jected to a fluid with a known velocity. Specially designed apparatus is normally 
used to generate high quality flow with uniform velocity, temperature and very 

low turbulence level. Hot-wire to be calibrated is kept firmly in the calibration 

apparatus such that the direction of velocity is normal to the wire. The flow 
should be steady. The voltage output of the wire is noted at zero velocity and at 
higher fluid velocities. The calibration data should be repeatable with very low 

scatter. 



6.1 Velocity Data Reduction 


102 


In the present study, calibration has been performed in the wind tunnel itself 
due to its high quality. The maximum turbulence level is 0.5% in the upstream 
section of the apparatus. 


Isothermal Calibration of Hot-wire 

DANTEC hot-wire equipment uses a linearizer card (56 N 21) which produces 
a voltage signal linearly proportional to flow velocity. The mathematical curve 
used by this card is given below. 


y = + D (6.1) 


Here A, B,C,D and E are constants. Normally constants are selected to give a 
fluid velocity equal to lOm/sec for a wire output of 10 V, Here y is the linearizer 
output equal to normalized velocity of fluid defined as 

y = 10 X (6.2) 

Umax 

and X is normalized voltage defined as, 


a: = 10 X 


JlziL' 

Vjnax ~ ^0_ 


(6.3) 


where Vo is the output of the wire measured at zero velocity of fluid and Vmax is 
the measured voltage at the maximum fluid velocity Umax- 

The constants A, B,C,D and E are obtained by means of an iterative least 
square error approach. For any calibration point (xi, yi) Equation 6.1 will produce 
an error (j)i given by 


+ (6.4) 

The least square approach requires the sum of errors at all calibration points to 
be minimum, i.e. 

4 = X = E ~ [10^'-®“'"®’“ + CXi + D] (6.5) 

0 . ( 6 . 6 ) 


and 


^ 


Let I be any of the parameters A,B, C,D or E, then 


iz=:l 


(6.7) 


Hence yl, S, C, D and B can be calculated using the above system of equation. 
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6.2 Temperature Data Reduction 

Temperature data was collected by RTD for which the resistance is determined by 
a precalibrated thermistor. The HP 3457A multimeter was used for the display 
of the resistance of probe wire. It was programed for converting the operating 
resistance of the probe wire into the corresponding temperature. MATH and 
STAT operations were utilized for the purpose of conversion. 


6,3 Liquid Crystal Thermography 

Thermocromic liquid crystals were spread uniformly on the test surface and after 
the rib. The liquid crystal experiment was performed in two phases. 

1. Calibration test to relate color with temperature. 

2. Transient test with two different Reynolds numbers to capture the colour 
images in a time sequence, around the rib. 


6.3.1 Calibration Test 

Prior to the unsteady temperature measurement over the test surface, the method 
of color expression was defined and the relationship between the color and tem- 
perature of the liquid crystal was assessed in a calibration test. Although the 
colors of the liquid crystals are observed in terms of color components RGB, the 
R, G and B-data has not been used to calibrate the sheet. Instead, the RGB color 
space is transformed to the hue-intensity-saturation space, because the hue value 
is a monotonic function of the dominant wavelength reflected by the ciyst js. It 
is therefore best suited for a unique color-to temperature relation. A perfectly 
lighted surface showed uniform background intensily of the color over the entire 
test surface. The background intensity determines the clarity of the picture seen 


Tte^trarface was heated to a temperature of nearly 43»C. 

of one thermocouple was monitored continuously by a PC “ 

j mi, liPfltpr DOwer supply was cut off. The calibration 

data acquisition card. Then the heater power supp , ,,,, to 

. a The plate Started to cool owing to 

tp<?f was done under the no flow condition. p 

test was done un ora+„rpd at every 0.2‘’C temperature interval. 

natural convection and images were captured at every 0 
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Images were captured until the color of the surface became red. Assuming that 
temperature is same within a small area, 45x45 pixels were selected in the image 
around the location where the thermocouple was placed. These images were 
processed in the PC to get HSI values from RGB values separately for those 
particular pixels. The HSI values were averaged to get a single hue, saturation 
and intensity value appropriate for the temperature. With appropriate curve 
fitting a plot was obtained between hue and temperature. Figure 6.1 shows the 
calibration curve which was derived from the transient experiment. It can be 
observed that intensity is almost constant for the entire range of temperature. 
Figures 6.3 and 6.4 show RGB and HSI plots in present calibration. 



Figure 6.1: Calibration chart of the liquid crystal used in the present work. 
Using hue from the HSf Spa« reduces possible sources of uncertainty, such 
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as the brightness of the light source [Farina et al, 1994] The hue is observed to 
have a monotonic change with temperature and also remains independent of the 
local illumination strength. The variance of hue with temperature is found to 
be small and constant, characterizing the uniform distribution of color over the 
chosen pixels. 


6.3.2 Transient Test 


The transient heat transfer method offers a significant advantage of yielding the 
local heat transfer coefficient over the complete test surface in a single experiment. 
The technique involves heating of the test plate to a certain temperature. The 
flow is then turned on. The lighting arrangement and the camera position are 
kept identical to that of the calibration test. After the plate reaches a particular 
temperature, the heater power supply is cut off while the flow continues to cool 
the surface. Due to forced convection the plate stars to cool. Images are then 
captured at 2 seconds time intervals. This experiment was carried out for two 
Reynolds numbers, namely 160 and 260. Here, the Reynolds number is based on 
rib height. Using the calibration curve, temperatures at different points over the 
surface were evaluated. 

During the transient heat transfer experiment each pixel location will reach a 
prescribed temperature at different times, depending upon the local heat transfer 
coefficient. The transient temperature distribution can be used for evaluation of 
the local heat transfer coefficient by assuming the test surface as a semi-infinite 
solid. For a semi-infinite solid, the one dimensional transient conduction equation 
is: 


dt 

Using a convective boundary condition at the surface, 

-k&\,,^ = hlT{v = 0,t)-U 


( 6 . 8 ) 


(6.9) 


where is the fcee-stream temperature (or the local bulk mean temperature 
(n„.) in the case of a duct) and T(t) is the temperature of the surface at any 
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time t. Assuming Too or Tbm to be constant, a semi-infinite wall with constant 
properties, and a uniform initial temperature of T^i, the solution for the temper- 
ature at the surface {y = 0 ) is 


or 


Tnji - T{t) 

Tyji Too 


= 1 - exp 




( 6 . 10 ) 


T{t) - T„ 
T'wi ~ Too 


= exp(7^)erfc(7) 


where the parameter 7 is defined as 

- 

\/pck 


( 6 . 11 ) 


( 6 . 12 ) 


The only unknown in the Equation 6.12 is h, which is determined by a least 
squares fit through the experimental data in a prescribed total time span At. 
It was observed that the changes in temperature with time at any particular 
pixel starts diverging from the theoretical profile based on the semi-infinite solid 
assumption. So, the total time span (At) for which the curve fitting is done, 
plays an important role. If At is very high, then at some pixels, where heat 
transfer rate is very high, would reach to a temperature which is beyond the color 
display range of the liquid crystal. At the same time, some other pixels where 
heat transfer is low, would reach to a condition where experimental temperature 
gradient starts diverging from theoretical curve. At these pixels the predicted 
heat transfer coefficient would be comparatively higher than the actual. Hence, 
the pixels with higher heat transfer rate give the estimate of the total time span 
(At), for which the curve fitting should be done at each pixel. 

Validity of Semi-iiifimte Solid Assumption 

The validity of the semi-infinite solid assumption is verified from the transient 
temperature measurement of the thermocouples located at the top and bottom 
of the bakelite plate and its thickness. The assumption of a semi-infimte solid is 
valid only if the specimen is made from a material with a low thermal diffusivity 
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Figure 6.2: Verification of semi-infinite solid assumption with temperature mea- 
surement at different locations of the plate 


and chosen to be sufficiently thick. Chan et al. [2001] suggested a criterion for 
the minimum thickness of the material {y), as: 

y > 4'N/aAt (6.13) 

In the present work, bakelite is assumed to be. the semi-infinite solid for which 
the thermal diffusivity is very low. The required thickness for the plate for even 
150 seconds of experimentation would be 16,54 mm. Since the thickness of the 
upper bakelite plate is more than 25 mm, the above condition (Equation 6.13) is 
fully satisfied. Temperature variation of the thermocouples placed at the top and 
bottom surface of the bakelite plate are shown in Figure 6.2. It can be observed 
that the temperature of the lower surface of the Bakelite is unchanged while the 
temperature of the top surface decreases with time. This observation confirms 

the validity of semi-infinite solid assumption. 

The other assumption in the present study is that conductivity of the alu- 
minum plate is very high and thickness is very smaU, so that the temperature 
gradient in the aluminum plate is equal to the temperature gradient of the upper 
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surface of the bakelite sheet. 

6.3.3 Inverse Approach for Nusselt Number Estimation 

An inverse forced convection formulation for transient heat transfer has been dis- 
cussed in Section 4.4. This approach has been adopted in the present work to 
obtain local Nusselt number from unsteady measurements of the plate tempera- 
ture. The calculation proceeds with the following assumptions: 

1. The velocity and the temperature fields are two dimensional. 

2. The flow field is steady. 

3. As the rib height to channel height ratio is quite small, being equal to 1/24, 
the free slip boundary condition has been applied in the far field of the rib. 
This assumption has been validated by suitable numerical experiments. 

4. Constancy of the local Nusselt number with time has been assumed. This 
is equivalent to requiring that the thermal field is fully evolved in time. In 
other words, the temperature in the fluid medium fall continuously with 
time during the cooling of the plate. The cooling rates are thus taken to be 
dependent only on the heat transfer coefficient. 

5. The effects of free convection and radiation are neglected. 
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Chapter 7 


Convective Heat Transfer from a 
Flat Plate 


This chapter presents the outcome of applying the inverse technique to predict 
the local Nusselt number for flow over a flat surface. Three different solvers^ are 
used namely, (1) steady two dimensional parabolic, (2) steady three dimensional 
parabolic and (3) unsteady three dimensional elliptic. Detailed formulations are 
furnished in Chapter 3 and the numerical issues are discussed in Appendices A and 
B. The goal of the simulation is to asses the capability of the inverse procedure to 
accommodate noise in the measured data and its ability to use additional (though 
limited) use high quality data. Study proceeds along the following steps: 

1. Solve direct problem assuming a meaningful plate temperature distribution. 
The plate temperature is steady but may vary with position over the plate. 
This step will generate interior fluid temperature field and the local Nusselt 
number variation over the surface. The solution of this problem is termed 
as exact solution. 

2. Select interior temperature on certain planes for inclusion in the functional. 

3. Perturb the plate temperature and the interior temperatures of step 2 with 
zero mean Gaussian noise (Appendix C). Thus each perturbed temperature 
fall in the range 

^(perturbed ~ ^lexact ^ C^-1) 

^The solvers refer to the numerical algorithm for solving the direct problem embedded in 
the inverse technique. 
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where a is the standard deviation of the Gaussian distribution. F indicates 
the confidence interval. For present analysis a 99% confidence interval is 
maintained by keeping -2.576 < F < 2.576. It is expected in practice that 
the standard deviation of the interior temperature is much smaller than 
that of the surface temperature. 

4. Filter the plate and interior point temperatures by Gram orthogonal poly- 
nomial smoothing. 

5. Solve the complete inverse problem to estimate the local Nusselt number 
distribution. 

6. Compare the Nusselt number distribution with the exact solution of step 1. 

Sensitivity analysis is presented in all cases to authenticate the behavior of the 
inverse algorithm. The sensitivity coefficient at any point of the domain defined as 
the absolute change in the temperature for a unit change in local Nusselt number. 
Values of the sensitivity coefficients are obtained from the numerical solution of 
the sensitivity equation. While the absolute values of sensitivity carries little 
physical interest, they refiect the relative influence of the wall boundary condition 
over the field data. 

. 7.1 Two Dimensional Steady Parabolic Solver 

The two dimensional parabolic problem is less complex and computationally less 
intensive. Hence it permits numerical experimentation with this formulation. 

The effectiveness of inverse approach in local Nusselt number estimation is 
shown in Figures 7.1, 7.2 and 7.3. The inverse solution has been compared with 
the direct solution given by filtered data as well as the exact solution. The quality 
of inverse estimation shown here is quite superior to simple filtering. 

In the present work it is observed that the success of any inverse algorithm 
depends upon sensitivity, quantity and the noise level in the input data. These 
factors are studied and discussed below. 
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7.1.1 Sensitivity Study 

It is to be realized that for the present problem the boundary information is 
contained only in the vicinity of the plate. Figure 7.4 shows the results of the 
sensitivity study. The maximum value of sensitivity occurs at the plate for a 
Reynolds number of 500, while the minimum value is zero. All the data has been 
normalized between these two values. The high sensitivity zone is observed to 
become smaller with an increase in the Reynolds number. This clearly conforms to 
the physical nature of the problem. The thickness of the zone which contains the 
boundary information, i.e. the thermal boundary layer thickness, becomes smaller 
with the increase of Reynolds number^. This is one of the major constraints in 
the design of an experiment in which additional data is recorded. 

As the measurement plane shifts from the surface, the quality of inverse esti- 
mation degrades. This is shown in Figure 7.5. This effect is studied by the noise 
free inputs from five different measurement planes. At each plane all tempera- 
tures are assumed to be known but the surface temperature information is not 
included. At higher Reynolds number it is observed that their is no improve- 
ment in the initial guess if the plane is placed beyond a critical distance. The 
sensitivity coefficients at this planes are found to be zero. These observations are 
self-consistent. 

7.1.2 Effect of Noise Level and Grid Fineness 

As the scatter in the data poses the primary constraint, any inverse analysis is 
considered incomplete if it does not account for it. In the present study we have 
considered two different error levels in data. The plate temperature is known 
at all points but with substantial errors, while a smaller amount of high quality 
interior temperatures are available. 

While very few noise-free interior data points can absorb a substantial noise in 
the plate temperature, a small noise in the interior data can degrade the solution 
to a large extent. These observations are presented in Figures 7.6 and 7.7. These 
trends confirm the ill-posedness of the inverse problem. During the study of 
scatter in plate temperature, it is assumed that 20% error-free data is available. 
The effect of interior data error has been studied when the surface data is not 


^Prandtl number is constant here 
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available. 

The negative effect of interior data error increases with Reynolds number. 
This fact can be explained from the sensitivity study (Figure 7.4). It shows 
that the sensitivity coeiRcient at any plane decreases with Reynolds number. For 
higher Reynolds numbers, the data quality deteriorates by both low sensitivity 
and scatter in the interior point temperature. 

The effect of data error are more critically examined by varying the amount 
and sensitivity of data. Figure 7.8 shows that diminishing the data quantity, data 
sensitivity and increasing data error, are all detrimental to the final solution. But 
their interdependence is quite subtle Figure 7.9 shows that a higher amount of 
data might not always be a guarantee of an improved solution, especially if the 
data error is substantially high. In fact, in some cases a lower amount of data can 
generate a better solution. With the increase of the number data points, noise as 
well as information are absorbed in the predicted solution. Which of these two 
will ultimately dominate is problem-specific. Higher scatter in input data from a 
higher sensitivity region may deteriorate the solution further. 

Unlike the direct problem, the solution of the inverse problem may not always 
tend to exact solution with the improvement in grid fineness, particularly if the 
data contains substantial scatter. This behavior is shown in Figure 7.9. Here the 
solutions have been derived on two different grids. In the finer grid, the number 
of additional data points is varied. It is observed that the solution On the finer 
grid with larger amount of data leads to a poor solution. Conversely the solution . 
quality improves with a reduced number of data. 

With the improvement in grid fineness, the change of the computed variables 
(i. e., temperature, in this case) decreases in between two consecutive grid points. 
If the standard deviation of data error becomes greater than the change in tem- 
perature in consecutive grid points, the inverse algorithm receives misinformation, 
causing degradation in accuracy. This phenomenon poses a great restriction on 
fix the optimum number of grid points for the inverse solution. 

7.1.3 Studies on Initial Guess and Uniqueness 

An effective optimization algorithm is expected to eliminate the remnants of 
initial guess from the final solution. However, this is not always true for inverse 



7.1 Two Dimensional Steady Parabolic Solver 


115 


problems, as the property of uniqueness may not be strictly applicable here. 
Figure 7.11 confirms this observation. If the experimental data exactly satisfies 
the energy equation, the final value of the functional would be zero; otherwise it 
will converge to a non zero positive value. If the amount of input data is very less, 
it can always be fit into the energy equation to produce an almost zero functional. 
If the data amount is large, the functional can be brought to a near zero value 
only when the data is exact. For a large amount of noisy data, irrespective of 
the boundary condition, it is to be expected that the solution will not conform 
point-wise to the energy equation. 

7.1.4 Studies on Stopping Criteria 

In the present study several stopping criteria have been employed. A comparative 
study is shown in Figure 7.12. If only one level of data scatter is assumed, all 
the criteria have shown similar performance. For two different scatter levels 
discrepancy principle is not applicable. It is realized that as the scatter level 
mismatch increases data smoothing produces better solution than the additional 
data criterion. 
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Figure 7.1: Comparison of inverse and direct solution with perturbed and filtered 
data. Scatter in plate temperature is .cr=0.2. Scatter in interior temperature is 
(7=0.0001, Only 50% interior temperatures are available at y = 0.01. Velocity 
and temperature have been obtained from a parabolic formulation. 
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Figure 7.2: Comparison of inverse and direct solution with perturbed and filtered 
data. Scatter in plate temperature is cr=0.2. Scatter in interior temperature is 
17=0.0001. Only 50% interior temperatures are available at y = 0.01. Velocity 
and temperature have been obtained from a parabolic formulation. 
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Figure 7.3: Comparison of inverse and direct solution with perturbed and filtered 
data. Scatter in plate temperature is cr=0.2. Scatter in interior temperature is 
(7=0.0001. Only 50% interior temperatures are available at y = 0.01. Velocity 
and temperature have been obtained from a parabolic formulation. 












Figure 7,4; Variation of sensitivity coefficient with distance from, the plate at 
different Reynolds numbers, (a) Variation of sensitivity at a given plane with 
Reynolds numbers. 
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Figure 7.5; Effect of the choice of the additional high quality data on inverse 
estimation. This data lies on the plane 2/= constant. 
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Figure 7.8: EflFect of the quality (standard deviation) , quantity and sensitivity of 
interior data on the inverse solution. 
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Figure 7.9: Effect of grid refinement on inverse estimation of the local Nusselt 
number. 
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Figure 7.12: Effect of stopping criteria on the inverse solution of local Nusselt 
number. 
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7.2 Three Dimensional Steady Parabolic Solver 

The effectiveness of the inverse formulation is revealed for three dimensional ge- 
ometries in Figures 7.13-7.16. The sensitivity analysis has also been carried out 
and presented in Figures 7.17 - 7.19. 

These figures clearly show that the inverse formulation presented here is ca- 
pable of generating meaningful results for the three dimensional thermal field. 
Specific observations are: 

1. Reasonable solutions have been achieved upto 10% scatter in interior data. 

2. For noise level and measurement plane remaining same, inverse estimation 
quality degrades with the increase in Reynolds number, 

3. The sensitivity analysis complies with the inverse estimations. 
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Figure 7.13: Comparison of inverse and direct solution. Scatter in plate temper- 
ature is (7=0.5, Scatter in interior temperature varies from cr=0.001 to c7=0.1 
. Only 50% interior temperatures are available at y = 0.01. Velocity and tem- 
perature have been obtained from a parabolic formulation. Velocity field is two 
dimensional. 
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7.3 Three dimensional unsteady elliptic formu- 
lation 

The performance of the inverse formulation is studied for three dimensional un- 
steady heat transfer from the flat plate (Figures 7.20 - 7.22). The intermediate 
direct problems have been numerically solved by an elliptic formulation^. The 
results from sensitivity analysis are shown in Figure 7.23. 

The inverse formulation is shown to be less effective for small time but predicts 
well as the time progresses. The sensitivity studies show that the sensitivity 
coefficient at a point increases with time. This allows the inverse algorithm to 
extract more and more informations from a point as the time progresses. The 
specific observations of this part of the study are; 

1. Unsteady inverse technique becomes more eflicient as time progresses. 

2. For noise level and measurement plane remaining same, inverse estimation 
quality degrades with the increase in Reynolds number. 

3. The sensitivity analysis complies with the inverse estimations. 


parabolic solver in space and time could also have been used, but would be inapplicable 
for flow past a rib that is discussed in chapter 8. 
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Figure 7.21; Comparison of inverse and direct solution. Scatter in plate temper- 
ature is cr=0.5. Scatter in interior temperature varies from cr=0.001 to <7=0.1 . 
Only 50% interior temperatures are available at y = 0.01. The temperature field 
has been obtained from an elliptic formulation. A two-dimensional parabolic 
solver is used for the flow field. 
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Chapter 8 

Convective Heat Transfer from a 
Flat Plate with a Rib 


This chapter presents the outcome of applying the inverse technique to predict 
the local Nusselt number for flow over a surface with a solid rib mounted over it. 
Three different solvers are used namely: (1) Unsteady two dimensional formula- 
tion with a steady plate temperature, (2) Unsteady three dimensional formula- 
tion with a steady plate temperature, (3) unsteady two dimensional formulation 
with experimentally determined time varying surface temperature. Detailed for- 
mulations are furnished in Chapter 4 and the numerical issues are discussed in 
Appendices A and B. Inverse solutions of (1) and (2) contain the following steps: • 

1. Solve direct problem assuming a meaningful plate temperature distribution. 
The plate temperature is steady but may vary with position over the plate. 
This step will generate interior fluid temperature field and the local Nusselt 
number variation over the surface. The solution of this problem is termed 
as exact solution. 

2. Select interior temperature on certain planes for inclusion in the functional. 

3. Perturb the plate temperature and the interior temperatures of step 2 with 
zero mean Gaussian noise (Appendix C). Thus each perturbed temperature 
fall in the range 

^(perturbed “ ^lexact ^ (®'l) 

where a is the standard deviation of the Gaussian distribution. P indicates 
the confidence interval. For present analysis a 99% confidence interval is 



8.1 Two Dimensional Formulation 


142 


maintained by keeping -2.576 < F < 2.576. It is expected in practice that 
the standard deviation of the interior temperature is much smaller than 
that of the surface temperature. 

4. Filter the plate and interior point temperatures by Gram orthogonal poly- 
nomial smoothing. 

5. Solve the complete inverse problem to estimate the local Nusselt n um ber 
distribution. 

6. Compare the Nusselt number distribution with the exact solution of step 1. 

For the transient plate temperature problem, the noise level of the experi- 
mental data is not known exactly. The stopping criterion based on additional 
data is used to stabilized the solution. No data smoothing technique has been 
incorporated for this part of the study. 

8.1 Two Dimensional Formulation 

The eflPectiveness of the inverse formulation is revealed for the two dimensional 
geometry in Figures 8.1 - 8.3, The results of the sensitivity analysis has also been 
presented in Figure 8.4. 

The sensitivity at the interior points is higher than for the flat plate case. 
This is attributed to the presence of rib, that influences the flow field to a greater 
extent, than a boundary-layer. This phenomenon serves as the reason for better 
prediction of the unsteady local Nusselt number as compared to a flat plate. 
Specific observations in this section are: 

1. Transient prediction of the solver is superior than the flat plate one. 

2. For noise level and measurement plane remaining same, inverse estimation 
quality degrades with the increase in Reynolds number. 

3. The sensitivity analysis complies with the inverse estimations. 
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Time Step=50 
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Figure 8.1: Comparison of inverse and direct solutions. Scatter in plate tempera- 
ture is a=0.5. Scatter in interior temperature is varying from .0001 to .01. Only 
50% interior temperatures are available at y — 0.01. Velocity and temperature 
have been obtained from an elliptic formulation. 
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Figure 8.2: Comparison of inverse and direct solutions. Scatter in plate tempera- 
ture is cr=0.5. Scatter in interior temperature is varying from .0001 to .01. Only 
50% interior temperatures are available at y = 0.01. Velocity and temperature 
have been obtained from an elliptic formulation. 
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Figure 8.3: Comparison of inverse and direct solutions. Scatter in plate tempera- 
ture is (j=0.5. Scatter in interior temperature is varying from .0001 to .01. Only 
50% interior temperatures are available at y = 0.01. Velocity and temperature 
have been obtained from an elliptic formulation. 
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8.2 Three Dimensional Formulation 

The effectiveness of the inverse formulation is revealed for the two dimensional 
geometry in Figures 8.5 - 8.7. The results of the sensitivity analysis has also been 
presented in Figure 8.8. 

The sensitivity at the interior points is higher than for the flat plate case. 
This is attributed to the presence of rib, that influences the flow field to a greater 
extent, than a boundary-layer. This phenomenon serves as the reason for better 
prediction of the unsteady local Nusselt number as compared to a flat plate. 
Specific observations in this section are: 

1. Transient prediction of the solver is superior than the flat plate case. Be- 
havior of the three dimensional solver matches with the two dimensional 
one. 

2. For noise level and measurement plane remaining same, inverse estimation 
quality degrades with the increase in Reynolds number. 

3. The sensitivity analysis complies with the inverse estimations. 



8.2 Three Dimensional Formulation 


148 


Exact Solution 




Inverse Estimation 





Figure 8.5: Comparison of inverse and direct solutions. Scatter in plate temper- 
ature is (7=0.5. Scatter in interior temperature is cr=0.01. Only 50% interior 
temperatures are available at y = 0.01. Velocity and temperature have been ob- 
tained from an elliptic formulation. The flow field is two dimensional with three 
dimensional temperature field. 
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Figure 8.7; Comparison of inverse and direct solutions. Scatter in plate temper- 
ature is (7=0.5. Scatter in interior temperature is <7=0.01. Only 50% interior 
temperatures are available at y = 0.01. Velocity and temperature have been ob- 
tained from an elliptic formulation. The flow field is two dimensional with three 
dimensional temperature field. 
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Figure 8.8: Variation in sensitivity coefficient with time and Reynolds number 
for inverse estimation of local Nusselt number over a ribbed surface. 
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Figure 8.9; LOT image at time 0, Re=260. 



Figure 8.11: LCT image after 45 sec, Re=260. 
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Figure 8.12: Local Nusseit number estimation using inverse technique with tran- 
sient LCT data. 
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Chapter 9 


Conclusions and Scope for Future 
Work 

9.1 Conclusion 

An inverse function estimation technique for two and three dimensional forced 
convective heat transfer from a flat plate has been discussed in the present work. . 
The focus of the research is towards determining the local Nusselt number distri- 
bution over the flat surface, and the effect of a surface mounted rib. The input to 
the calculation is the surface temperature variation as well as selected tempera- 
tures in the fluid phase. The study supports the conjugate gradient method as an 
efficient approach for the solution of inverse convection equations. The procedure 
also suppresses within limits the effect of initial guess in the final solution and 
scatter in the input data. Sensitivity study, required during the inverse procedure 
is seen to be a useful tool for designing an experiment. 

The inverse forced convection modeling has also been found to be useful for 
extracting information from noisy experimental data. The results show that 
the tool is promising and deserves further exploration. Certain advantages and 
disadvantages are to be noted regarding the conjugate gradient approach: 

Advantages 

1. The procedure is accurate and the exact solution is realized if the data given 
is exact. 

2. The procedure yields a qualitatively good prediction when noise in the 
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interior data is less than 0.1, namely 10%. 

3. The technique subdues substantial surface temperature scatter (0.4) pro- 
vided the interior data is of superior quality. 

4. The number of iterations needs to minimize the functional is not high, 
mostly falling in the range of 3 to 50. 

5. The procedure does not need a priori information about the unknown func- 
tion, the Nusselt number in the present context. 

Disadvantages 

The computational cost increases rapidly with the increase of problem complexity, 
Reynolds number, problem dimensions, and grid fineness. All computer programs 
were executed in a Pentium4 PC with 512 SDRAM. A comparison of CPU times 
is given in the following table (Reynolds number is 10000 for fiat plate and 260 
for the ribbed surface). 


Formulation 

CPU time 

Two dimensional parabolic, flat plate 

30 seconds 

Three dimensional parabolic, flat plate 

3 hours 

Three dimensional elliptic, flat plate (upto 50 time steps) 

2 days 

Two dimensional elliptic, ribbed plate (upto 50 time steps) 

4 days 

Three dimensional elliptic, ribbed plate (upto 50 time steps) 

9 days 

Two dimensional elliptic, transient plate temperature, 
ribbed surface (upto 19000 time steps) 

18 days 


The CPU time needed for the velocity computations are not included in the 
above table. The CPU time for velocity computation varies from 60 seconds (for 
two dimensional parabolic problem) to 4 days (vorticity-stream function formu- 
lation). 

9.2 Scope for Future Work 

The present work has opened up new areas of exploration. Following challenges 
are recommended for future researchers. 

1. An extensive assessment of the applicability of inverse forced convection 
techniques against life experiments. 
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2. The relative performance study of different inverse convection techniques 
in extracting useful informations from noisy experimental data. 

3. Acceleration techniques to decrease the CPU time. 

4. An inverse procedure for obtaining the skin friction coeificient from velocity 
data. 

5. Identification of dynamic filters for transient problems. 
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Appendix A 

Discretization and Numerical 
Algorithm 


A.l Solution of Velocity Field 


In the present work steady, laminar two-dimensional flow of incompressible, New- 
tonian fluid is addressed. Two different solution methodologies have been applied 
for computing the velocity fleld. 


1. Two dimensional boundary-layer approximation. 


2. Vorticity - stream function formulation. 

The appropriate differential equations have been solved by the method of 
finite diflferences. The numerical solution has been developed on uniform and 
non-uniform grids. 

A. 1.1 Two Dimensional Boundary Layer Flow 


Recalling the formulation presented in Chapter 3, we have the following governing 
equations for a fiat plate boundary-layer: 


d{u^) 

dx 

The boundary conditions are: 


du dv 
dx'^ dy 

(A.1) 

d{uv) 1 d'^u 

dy Re dy^ 

(A.2) 
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y 



y = 0, u — O, V = 0 
y — d, u = 1 


(A.4) 

(A.5) 


The momentum equation is solved by the Crank-Nicolson scheme, using central- 
difference discretization for both advection and diffusion terms. Hence 


- {u% 1 

Arc 2 
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2Ay 

i+1 _ 97/^+! 4- 7/^+^ 
■j+i 


2Re 
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'i+i 


2Ay 


A2/2 


(A.6) 


This leads to a tridiagonal system of equations: 

Au^J,+Bv^^^ + Cuy^\ = D (A.7) 

with the following boundary conditions (Figure A.l) 


y = Q, 

y = d, = 1 


The coefficients in Equation A.7 are 
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j;*+l 

4Ay 

i+l 
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Ax 
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1 

2ReAy^ 
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ReAy^ 


(A.8) 

(A.9) 


(A.10) 


(A.11) 
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y*+i 

C = ^ 


4Ay 2ReA2/2 


Aa: 4 


{nvYj+i - 

Ay 
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1 fui+^-2ui+u)_^ 


2Re 


Ay2 


(A.12) 


(A.13) 


Prom the continuity equation, we can determine the y component of velocity 


as 


v = -^Judy 


(A.14) 


In finite difference form, 

-3 // udy . + 4 udyj ^ - // udy _ 





2Ax 


(A.15) 


For 2 = 1 


and 


„i+l — 


Aa; 


J, = 0, tif ' = 0 


(A.16) 

(A.17) 


Thus knowing the values of u and v at the i-th position, we can find the values 
of u and v at the i + 1-th position. However, the non-linearity and the coupling 
of the above equations leads to an iterative solution comprising of the following 
steps: 


1. Set 2 = 1 and assign tt* and u* for all j from the given initial condition. 

2. Compute D. 

3. Assign iteration counter k = 0 and set[u*'''^]*^ = u*, = u*. 

4. Set and for all j to compute A,B,C 

5. Solve the tridiagonal system of equations using Thomas algorithm to find 

for all j. 

6. Find for all j using Equation A.13 or A.14. 

7. Set and [ 2 ;'+^]"+' = 2 ;'+' for all j. 
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8. Compute the global convergence quantity J 



9. If J < e (e is a small real number close to zero), 
set and 

10. Else, set it = A: + 1 and go to 4. 

11. Set i = « + 1 and go to 2. 

A. 1.2 Vorticity - Stream Function Formulation 


\l~*\ h~ ~H 


H 


^2 

In frl 


Figure A.2: Plate with Single Rib 


Recalling the formulation presented in Chapter 4, we have the following gov- 
erning equations for steady, two dimensional, incompressible flow: 


duj d{uu) d{voj) _ 1 / d^u) 

~d^ ~W “ Re 

dx^ dy^ 


(A.18) 

(A.19) 
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u 


V = 


(A.20) 


dip 
dy 

For a surface mounted rib (Figure A.2) the initial and boundary conditions are: 


t = 0, a; = 0, Ip = y, ^ = 1, V = 0 


(A.22) 


a: = 0, oj = 0, Ip = y, u=l, v = 0 
doo d'^ip du dv 


y = 0, u = —'^‘^ip, V’ = 0) w = 0, V = 0 
y = d, u = -V^ip, ip = y, u = l 
X = ^ 1 , and y < dn, u = —V^ip, ip = 0, u = 0, u = 0 
X = li + Iri, and y < dri, w = —V^ip, ip = 0, u = 0, v — 0 
X < li + Iri, and X > li, and y = dri, w = —V^ip, .'0 = 0, u = 0, u = 0 (A. 29) 
For nodes falling within the rib 

X < li+ Ifi, and X > li, and y < dri, u = 0, v = 0 (A.30) 


(A.23) 

(A.24) 

(A.25) 

(A.26) 

(A.27) 

(A.28) 


An explicit time-marching scheme is applied for the vorticity transport equa- 
tion. The Poisson equation for the stream function is solved in each time step 
using Gauss-Seidel iterative method with successive over-relaxation. The convec- 
tive terms are discretized by QUICK scheme (Quadratic Upstream Interpolation 
for Convective Kinematics) [Leonard, 2002, Minkowycz et al, 1988] and central 
differencing has been used for the diffusion terms. 

The discretized form of the vorticity transport equation can be expressed as 
follows. The following notation has been employed; superscripts and subscripts 
denote time and space co-ordinates respectively 


w, 


n-f 1 
hj 


<. + At 


Re I 


-At 


U, , I .OJ. , I , — Vj. 1 -U/. 1 . 
2+2^ I 2>J * 2^^ 

Ao: 


u 


+ 


ij+l ij+5 


Ay 


(A.31) 
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The boundary conditions for the surface mounted rib geometry are 

X ■ 


= 0, cuf/i = 0 


y. _ I ^n+1 _ 

*' ’ (As,)- 


s/ = i, a,”/* = 

a; = ^ 1 , and y < dri, oj 


i,i (Ay)2 

- 2s/,.y_\ + c-?j) 

(Aa;)2 


,/,n+l 


(Ax)^ 

X = + Zri, and t/ < dri, 

: < + Zri, and X > Zi, and y = dri, = ~ 

tream function equation can be discretized as 

i+ij ^ n,j-i _ , .w+i 

(ASP (Ay)2 « 

oisson equation is solved within 
opriate boundary conditions axe 


(A.32) 

(A.33) 

(A.34) 

(A.35) 

(A.36) 

(A.37) 

(A.38) 


(A.39) 

V — ^ / 

1 solved within every time step for the stream function, 
conditions are: 


For uniform incoming flow x = 0, = y 

(A.40) 

X = = 2v-rJ-i- - V'Ki 

(A.41) 

y = 0, = 0 

(A.42) 

y = d, = dy 

(A.43) 

X = Zi, and y < cZri, = 0 

(A.44) 

X = Zi + Zri, and y < cZ^i, = 0 

(A.45) 

X <li + lri, and X > Zi, and y = dri, tpij^ = 0 

(A.46) 


..n+1 _ 


(A.47)', 



A.2 Solution of Energy Equation 


171 


„n+l _ n-l.j Wi+\,j 


h3 


2i\.x 


x = o, u'yf = 1 , = 0 


X — Ip, u'lj 




■i-l,j n+1 _ ~ ^i-2,j 

J 


n+1 


,n+l 


y = 0, = 0, 

y = d, ug-^ = 1 

X < + lj. 1 , and X > k, and y < dri, w"7^ = 0, = 0 


(A.48) 

(A.49) 

(A.50) 

(A.51) 

(A.52) 

(A.53) 


The sequence of calculations for obtaining the velocity distribution has the fol- 
lowing steps: 


1. Specify initial values for lj and ^ at time t = 0 

2. Find u and u at t = 0 

3. Solve the vorticity transport equation for u) at all interior points at the next' 
time step. 

4. Solve Poisson equation for stream function within the time step. 

5. Find the velocity components u = V'y and v = —^ix- 

6. Determine the values of u at the boundaries. 

7. Return to step 2 for the next time step. 

A.2 Solution of Energy Equation 


In the present work, forced convection alone from the heated wall to the fluid has 
been considered. Three different strategies have been applied for the solution of 
the energy equation. 

1. Steady solution with two dimensional boundary-layer approximation. 

2. Steady solution with three-dimensional boundary layer approximation. 

3. Unsteady solution of the elliptic form of the energy equation. 
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y 



All the three methodologies have been applied to the flat plate geometry; 
only the third one is applicable for the ribbed plate. Finite difference method 
of discretization has been used to solve the differential equations. Both uniform 
and non-uniform grids are used. 


A.2.1 Two Dimensional Boundary Layer Solution 

Recalling the formulation presented in Chapter 3, we have the following governing 
equation for heat transfer from a heated flat plate: 


d{uff) d{vd) 1 5^0 

dx * dy Pe 

(A.54) 

The boundary conditions are: 


11 

O 

11 

O 

(A.55) 


(A.56) 

o 

II 

11 

(A.57) 

where Nu is the Local Nusselt number. 

The equation is solved by the Crank-Nicolson scheme, using central differences 
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for both advection and diffusion terms. Thus 


(««)■+* - M)* 1 


Ax 


+ 2 


- m-i 

2Ay 2 Ay 


1 ( 6f^\ - 26f^ + ^ 26] + ei _, ' 


2Pe 1 Ay2 

This leads to a triadiagonal system of equation 


Ay2 


Adit\ + Bep-^ + Cdp.\ = D 


with the following boundary conditions (A.3) 

y^O, = 1 




The coefficients of Equation A.59 are 


A = 




c = 


4Ay 2PeAy^ 
1 

PeAy2 
1 


— j 

B = — — ^ + 


Ax 

,.i+l 


4Ay 2PeAy2 


D = 


Ax 


Ay 


+ 


1 + 


2Pe 


A2/2 


(A.58) 


(A.59) 

(A.60)- 

(A.61) 

(A.62) 

(A.63) 

(A.64) 

(A.65) 


The above tridiagonal system of equation is solved by Thomas algorithm. 
Thus knowing the i-th values of 0, the i + 1-th values of 9 can be computed. The 
complete algorithm for the temperature field computations is as follows: 

1. Set i = 1 and assign 0* for all j from the given initial condition. 

2. Compute coefficients A, 5,(7, H. 

3. Solve the tridiagonal system of equations (A.59-A.60) using Thomas algo- 
rithm to find for all j. 


4. Set i = i + l and go to 2. 
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A. 2. 2 Three Dimensional Boundary Layer Solution 


Recalling the formulation presented in Chapter 3, we have the following governing 
equation for steady three dimensional flat plate boundary layer: 


dx 


+ 


The boundary conditions are: 


div9) 1 /d^e d^e\ 
dy ~Pe[dy^^d^) 

(AM) 

II 

0 

II 

0 

(A.67) 

n 96 

z/ = 0, — = -Nu 
dy 

(A.68) 

y = d, 6 = 0 

(A.69) 

0 

li 

0 

II 

(A.70) 

69 

(A.71) 


where Nu is the Local Nusselt number. In the inverse calculation, the local values 
of Nusselt number are updated iteratively. 

The equation is solved by the ADI scheme. The discretized form of the gov- 
erning equation is as follows: 


For i to i -I- 1: 


^ 2At/ 


1 

Pe 


For i 4- 1 to z -I- 1: 




(A.72) 


^ 2Ay 


1 

Pe 


- 2^^ + Ct , 


+ 


(A.73) 


(Ayf ' {Azf 

Equation A.72 produces tridiagonal system of equation in the y direction. For 
the frame at z -b 5 we can have the solution by sweeping in the z direction and 




Figure A.4: Heat Transfer over a Flat Plate 

solving tridiagonal matrix at each 2 location. Similarly, Equation A. 73 produces 
a tridiagonal system of equation in the z direction. So for the frame at i + 1 we 
can have the solution by sweeping in y direction and solving tridiagonal system 
of equation at each y location. 

A. 2. 3 Unsteady Energy Equation Solution 

Recall the formulation presented in Chapter 3, we have the following governing 
equation for unsteady three dimensional flow over a heated surface (Figure A.2, 
A.4 ). 


dO diuO) d{vO) 1 dH dH\ 

dt'^ dx dy “ Pe dy^ dz^ ) 

The initial and boundary conditions for the geometry of Figure A.2 are: 

(A.74) 

0 

II 

0 

11 



(A.75) 

a: = 0, = 0 

ox 



(A.76) 

CC = 1, 0 = 0 



(A.77) 

. dd __ 

y = 0, -^ = -Nu 
dy 



(A.78) 

0 

II 

^ i 

II 



(A.79) 

z = 0, ^ = 0 



(A.80) 
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2 = r, = 0 (A.81) 

where Nu represents the local Nusselt number distribution over the heated bound- 
ary. 


A fully implicit finite difference procedure has been employed. The convec- 
tive terms are discretized bu the QUICK[33,29] algorithm. Central differencing 
has been applied for the diffusion terms. The discretized form of the governing 
equation is as follows: 


7/"+! _ 7/” 

dt 


+ 




n-hl 


vn-f 1 


Aa; 


+ 


Ay 


^z+l,j,fc ^%j,k + ’^i-l,j,k \j+l,k ^%},k + \j-l,k ^i,j,k+i ~ ^%j,k + 

PeAx^ PeAy2 PeA^^ 

(A.82) 


The initial and boundary conditions for the geometry of the figure can be written 
as 


o 

II 

o 

II 

(A.83) 

11 

o 

II 

o 

(A.84) 

^ = 1. + Kli = 0 

(A.85) 

0. KJ:i - = 2A!,n<« 

(A.86) 

V = d, = 0 

(A.87) 

z = o, = o 

(A.88) 

^ = >-. = o 

(A.89) 


Equations have been solved by G-S iteration. The solution methodology consists 
of the following steps [Roche, 1976, Tannehill et al, 1997]: 

1. Specify initial values for at time t = 0 

2. Solve equations A.79 to A.86 for at time t = t + At using the Gauss- 
Seidel iterative method. 

3. Continue marching in time till t = t-hAt reaches its specified limit or the 
change in the local two successive time steps is negligible. 
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A.3 Solution of the Adjoint Problem 

The adjoint problem differs from the direct in the following sense. 
1. It contains a source term. 


2. It contains a negative diffusion term. 

3. It contains a final condition instead of an initial one. 


Let us consider a typical adjoint equation arising from the three-dimensional 
unsteady temperature field. Recalling the formulation presented in Chapter 3, 
we have the following adjoint equation and boundary conditions: 


t = ff, A = 0 
x = 0, A = 0 
X = 1, ^ + uXPe = 0 

2 / = 0, ^ = —2wPe{9u 
ay 

y = d, A = 0 


Tu;) 


dz 


0 


(A.91) 

(A.92) 

(A.93) 

(A.94) 

(A.95) 

(A.96) 

(A.97) 


Letting r = tf - 1 we have the following transformed system: 


dr dx dy 


'd^x d^X d^X 

Pe V dx^ dy^ dz"^ 


r = 0, A = 0 


x = l, 


dx 


+ 2{e-Y)5{y-yx) (A.98) 

(A.99) 
(A.IOO) 

+ uAPe = 0 (A.lOl) 


X = 0, A = 0 
dX 


y = o, 


dy 

y 


: —2wPe{6^ — Y^) 
d, A = 0 


(A.102) 

(A.103) 
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n 

dz~^ 

(A. 104) 

dX 

(A.105) 


The equation is now reduced to an initial value problem. The solution method- 
■ ology of the direct problem can now be used. However, the negative advection 
terms and the Dirac Delta function must have to be treated carefully. 


A.4 Discretization of Convective Terms 


In the parabolic problems, the convective terms have been discretized by a central 
difference scheme. For the elliptic equations, two and three dimensional QUICK 
schemes have been used. The QUICK algorithm (Quadratic Upstream Interpo- 
lation for Convective Kinematics) is based on local upstream-weighted quadratic 
interpolation in estimating numerical cell face values and gradients of all field 
variables undergoing simultaneous convection and diffusion. For pure convection, 
the scheme is third order accurate. 

To elaborate the scheme, let us consider one dimensional convection and dif- 
fusion of a scalar its transport is governed by the equation 


dt ^ dx dx^ 


(A.106) 


where ti is the convective velocity in the x direction and F is the diffusion coeffi- 
cient. Considering a uniform grid the term ^ can be written as 



d(t> 

dx Ax 

(A. 107) 

now 

= '^IN - icURVAx^ 

(A.108) 

where 

<^LIN = 2 

(A.109) 

and 


(A.llO) 


cuRv = '^«'f 

i\x^ 2 

(A.lll) 
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For two and three dimensions (i.e. if the scalar is a function of all three dimen- 
sions) the Equation A. 105 becomes 


= ^LiN 

- + 1 

CURVTYAy2 

(A.112) 

and 




•i’i-W = '^LIN - gCURVAi^ + icURVTYAs 

^ + A CURVTZ A.^2 

24 

(A.113) 

where 




CURVTY = 

— 2^j_i j -f 

Ay2 

if Ui_y > 0 

(A. 114) 

CURVTY = ^ 

i,j+l ~ 

Ay2 

o 

V 

1 

(A.115) 

CURVTZ = “ 

Az2 

, if > 0 

(A.116) 

CURVTZ = ' 

Az2 

O 

V 

1 

(A.117) 


This discretization scheme of the convective terms has been employed in all 
elliptic direct, adjoint and sensitivity problems and in the vorticity transport 
equation. The scheme is found to be free of any oscillation for the range of grid 
sizes, Reynolds numbers and Peclet numbers studied. 


A. 5 Time Step Estimation 


For accuracy of numerical simulation, the mesh size must be chosen small enough 
to resolve the expected spatial variations in all dependent variables. Once a mesh 
has been chosen, the condition on the time step necessary to prevent numerical 
instabilities are determined from the combination of Courant-Friedrichs-Lewy 
(CFL) criterion and the restriction of the grid Fourier number. According to the 
CFL condition, the distance the fluid travels in one time increment must be less 
than one space increment. This leads to a constraint on the time step in the form 


At < 


Ax Ay 


(A.118) 


The grid Fourier number restriction can be applied as 


At < 


{AxAyy 


Re 


2 [(Aa:)2 + (Ay)^ 


(A.119) 
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The minimum of this two is the upper limit of time step for explicit compu- 
tation. For implicit formulation the time step is taken less than 10 times of the 
corresponding explicit time step. 


Appendix B 

Grid Independence and Code 
Validation 

B.l Grid Independence Test 

The performance of any numerical scheme with respect to accuracy is highly 
dependent on the size of grid chosen. In the present study extensive numerical 
experiments have been performed to finalize the grid size for computation. It is 
observed that as the Reynolds number increases, the grid size has to be made 
smaller to have a grid independent solution. 

For the flat plate, the domain size is 1 x 1. The Reynolds number range is 
500 to 100000. A uniform grid is used with mesh size varies from from .01 x .01 
to .001 X .001 has been used. 

For ribbed surface, the domain length is 47.5 x 6.5. The grid is non-uniform, 
varying in a geometric progression. For a Reynolds number up to 845, it is 
observed that a grid points of 251 x 71 is sufficient. The same fineness has been 
maintained in the present work for all Reynolds number, ranging from 100 to 
300. The grid size varies from .01 close to the rib to .6 away from it, as shown in 
Figure B.l. 

A grid independence test result for the two dimensional energy equation is 
shown in Figure B.2. The Reynolds number for this plot is 300 (based on the rib 
height). The local Nusselt number distribution has been plotted and three grid 
sizes are tested. 



Local Nusseit Number 
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Figure B.l: Grids for a Ribbed Surface. 



Figure B. 2; Grid independence test for two dimensional energy equation. 
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b. 2 Code Validation 

All the computer programs have been extensively validated against analytical, 
experimental and published numerical results. 

B.2.1 Two Dimensional Boundary Layer Solution 

Two dimensional velocity and temperature field can be solved analytically with 
the boundary-layer approximation. This gives the following expressions for skin 
friction coefficient and the local Nusselt number; 

C/ = 0.644Re;^ (B.l) 

Nu^ = 0.332Pr5Re| (B.2) 

The numerical solutions were compared with these analytical solutions and are 
shown shown in figure B.3 

B.2. 2 Three Dimensional Boundary Layer Solution 

To validate the three dimensional boundary layer solution of the energy equation, 
a constant temperature over the whole plate was utilized. The solution along the 
mid-plane (along direction) has been compared with the analytical and two 
dimensional boundary-layer solutions in Figure B.3 

B.2. 3 Vorticity-Stream Function Solution 

The two dimensional velocity field has been solved by the vorticity-stream func- 
tion formulation in a non-uniform grid shown of Figure B.l. The velocity code is 
validated against experimental results and published numerical simulations. 

The comparison is against the numerical solution of backward facing step 
problem [Gartling, 1990]. The Reynolds number of flow being based on channel 
height is 800. The backward facing step height is half the channel height. The 
inlet velocity profile is parabolic. 

The experimental validation has been conducted in the wind tunnel described 
in Chapter 5. Hot wire anemometry (single wire probe) is used to measure the 
streamwise velocity distribution. As the hot wire anemometer is unable to capture 
the direction of velocity, near rib region region shows a mismatches. 

The comparisons are presented in Figure B. 4. 
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B.2.4 Two Dimensional Elliptic Energy Equation Solu- 
tion 

The two dimensional elliptic energy equation solver has been compared with ex- 
perimental and published numerical results in Figure B.5. The numerical result 
has been compared with a laminar flow solution past a heat generating obstacle 
[Leung et al, 2000] at three different Reynolds numbers. Here the obstacle is 
placed in an insulated channel. The channel height is considered as the char- 
acteristic length. The zero gradient condition has been applied on the outflow 
plane. 

Experimental validation has been carried out in the wind tunnel described in 
Chapter 5. For this validation, flow is established over the plate while the heating 
is on. After a long time the temperature field becomes steady. Temperature 
profiles at selected locations have been measured using RTD, and compared with 
the numerical simulation [Figure B.5]. 

B.2. 5 Three Dimensional Elliptic Energy Equation Solu- 
tion 

Three dimensional elliptic energy equation is solved for constant plate tempera- 
ture and mid plane results are conxpared with two-dimensional energy equation 
solver: The results are shown in Figure B.6. 



Skin Friction Coefficient 
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Distance Distance 


Figure B.3: Comparison of analytical and numerical solutions of two dimensional 

momentum and two/three dimensional energy equation. 
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Numerical Comparison with Backward Facing Step Problem 

Re=800 

$ $ $ Present Solution 



Comparison with Experimental Results 
• • « « Experimntal 








Figure B 4: Validation of the stream function- vorticity solution. 
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Present Solution ^ ^ Published Solution — 
Comparison with Numerical Results 



Comparison with Experimental Results 








Figure B. 5: Validation of the two-dimensional energy equation. 
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Two dimensional soiution 

Three dimensionai soiution « « « « 








Figure B.6: Comparison of two and three dimensional elliptic energy equation 
solution. 










Appendix C 

Gaussian Random Number 
Generation 


The scatter in temperature measurement has been assumed to have a Gaussian 
probability distribution. For simulating experimentally derived scatter, the Box- 
Muller method is used to transform white noise to Gaussian noise [Press et al, 
2000]. White noise has a uniform variate and is obtained from built-in function 
RAND of the FORTRAN library. A brief outline of the Box-Muller method is 
given here. 

If a:i,a:2, . . . are random deviates with a joint probability distribution of 
p{xi,X2, . ■ .)dxidx2 . . . and if yi,y2 , ... are each function of all the xi,X2,... (same 
number of yi, j/2, ... bs Xi,X2, ■ ■ . then the Joint probability distribution of the 


yi,y2, ...IS 

p{yi, y2 , . ■ .)dyidy2 . • . = pixi, X2 , . . .) 


d(Xi,X2,. 


% 1 , 2 / 2 ,. 



dyidy2 . . . 


(C.l) 


where |5()/5()| is the Jacobian determinant of the xi, X2,... with respect to the 
Vi ) 2/2? ■ ■ • (or reciprocal of the Jacobian determinant of the yi, 2/2, • • . with respect 

to the a:ri,rr2) • • .). 

An important example of the use of equation (C.l) is the Box-Muller method 
for generating random deviates with a normal (Gaussian) distribution, namely 

p{y)dy- 




e 2 dy 


(C.2) 


Considering the transformation between two uniform deviates on (0,1), xi, X2 and 
two quantities yi, 1/2 we get 


= A/-21na;iCos27r2;2 


(C.3) 
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2/2 = \/— 2 In Xi sin 2 nx 2 

(C.4) 

Equivalently we can write 


xi = exp[~{yl + y^)] 

(C.5) 

1 2/2 

X 2 = — arctan — 

27r 2/1 

(C.6) 

Now the Jacobian can readily be calculated 


d{xi,x 2 ) ^ p p _ [1 yJ] r 1 41 

(C.7) 


Since this is the product of a function of y 2 alone and a function of yi alone, we 
see that each y is independently distributed according to the normal distribution 
(C.7). 

There is an easy way of implementation of this algorithm as shown by Press 
et al. [2000]. This calls for the following steps: 

1. Choose two uniform variate (ui- and V 2 ) within the unit circle (not unit 
square), so that their sum of squares = Vi+ u|) is within the circle. 

2. Define a uniform variate Xi = R^. 

3. Define Vi/v 2 = 27rx2 which gives cosine and sine in Equations C.3 and C.4 
as Vi/R and V 2 /R respectively. 

V ' ' : ■ 

4. Using Equations C.3. and C.4 find the corresponding Gaussian random num- 
ber ?/i and 2 / 2 - 




